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Abstract 

We aim to constrain the evolution of AGN as a function of obscuration using an X-ray selected 
sample of ~ 2000 AGN from a multi-tiered survey including the CDFS, AEGIS-XD, COSMOS and 
XMM-XXL fields. The spectra of individual X-ray sources are analysed using a Bayesian method¬ 
ology with a physically realistic model to infer the posterior distribution of the hydrogen column 
density and intrinsic X-ray luminosity. We develop a novel non-parametric method which allows us 
to robustly infer the distribution of the AGN population in X-ray luminosity, redshift and obscuring 
column density, relying only on minimal smoothness assumptions. Our analysis properly incorporates 
uncertainties from low count spectra, photometric redshift measurements, association incompleteness 
and the limited sample size. We find that obscured AGN with Nh > 10^^ cm“^ account for 771*15% of 
the number density and luminosity density of the accretion SMBH population with Lx > 10'^'^ erg/s, 
averaged over cosmic time. Compton-thick AGN account for approximately half the number and 
luminosity density of the obscured population, and 381*17% of the total. We also find evidence that 
the evolution is obscuration-dependent, with the strongest evolution around Nh ~ 10^*^ cm“^. We 
highlight this by measuring the obscured fraction in Compton-thin AGN, which increases towards 
z ~ 3, where it is 25% higher than the local value. In contrast the fraction of Compton-thick AGN 
is consistent with being constant at « 35%, independent of redshift and accretion luminosity. We 
discuss our findings in the context of existing models and conclude that the observed evolution is to 
first order a side-effect of anti-hierarchical growth. 

Subject headings: galaxies: active — surveys — quasars: supermassive black holes — X-rays: galaxies 


1. INTRODUCTION 

Supermassive Black Holes (SMBH) are abundant in 
the local universe " every nearby massive galaxy har - 
bours one (Richstone et al. |1998[[Kormendy fc Ho(2013[ ). 


The maiority of the growth of these SIVIBH must have 

^ - 


been through a ccretion processes (Soltan||1982 Merloni 
& Heinz 20081. The material accreted by the SMBH 


likely originates from galactic scales and therefore one 
might expect some relation between the properties of 
galaxies and their central black holes. This idea was re¬ 
inforced by the discovery of tight relationships between 


ies ( 

Magorrian et al.| 

1998 Ferrarese & Merritt 2000 

Gebl 

tiardt et al.||2U(JUl I'i'remaine et al.||2(JU2| |h'errarese| 

(k i^brd 20051 K.ormendy & Do 

2013p, which can be in- 


host galaxies. For instance, the growth of AGN through 
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accretion may be linked to the host galaxies star for¬ 
mation, either weakly through a shared gas reservoir, or 
strongly via so-called “feedback” processes. In the latter, 
the energetic output from the AGN into the galactic en¬ 


vironment can strongly distnr b the galaxy (|Silk & Rees 
1998 Fabian fc Iwasat^|1999 1. 


To nnderstand the growth of SMBH over cosmic time 
and the relationship of that growth to galaxy evolution, 
it is first necessary to estimate the space density of ac¬ 
creting SMBH, and its evolntion over cosmic time in an 
unbiased way. This is challenging becanse it is known 
that many AGN are enshronded in gas and dust, making 
them difhcnlt to detect directly. X-ray emission is an ef¬ 
ficient way to reveal these obscured SMBH, at least for 
moderate column densities. On the other hand when the 
column of equivalent neutral hydrogen exceeds unit op¬ 
tical depth corresponding to the Thomson cross section 
(^Nh ~ 1.5 X 10^^ cm“^) AGN become difficult to find 
even at X-ray wavelengths. In this so-called “Gompton- 
thick” regime the AGN can none the less still be identified 
by hard X-ray emission which can emerge from the thick 
covering and/or via radiation that is reflected and scat¬ 
tered into the line of sight. These processes result in a 
characteristic shape to the X-ray spectrum, including a 
flat continuum and intense iron Ka emission lines which 
can be identified nsing spectral analysis. The contribu¬ 
tion of these most heavily obscured AGN represents the 
major remaining uncertainty in our knowledge of the ac¬ 
cretion history, and placing limits on their contribution 
is therefore vital. 

Previons work has attempted to constrain the num- 
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her of Compton-thick AGN using X-ray background syn- 
Ddr 


thesis (e.g. Gilli et al.| 20071, but while this popula¬ 


tion is needed to reproduce the shape of the XRB spec¬ 
trum, the method is relatively in sensitive to the pre - 
cise Compton-thick AGN fraction (Akylas et al. 20121. 
Multi-wavelength data have also been exploited to iden- 
tify Compton-thick A GN and constrain their number. 
These include optical dRisaliti et al. 1999 Gappi et al. 


2006 JPanessa et al.| 2006 


2009 iGilli et al.pOlOjlVignali et al.|2010 I^Tia et al. 


Migiioli et al.||2bl3| I 


jAkylas fe Ceorgantopoulos 

a,li et al.|2010t|Jia et al.pOl^ 

ed (e.g’ 


ignali et al. 20141, infrared (e.g. 

[Alexander et al.|2011t|Brightman 

and hard X-ray (e.g. [Sazonov et 


& JNandra||2011a|b _ ___ 

20081 Burlon et ^[2011[ [Alexander et al.| 20R{HLanzuisi 


et al. 2014p diagnostics applied to local and non-local 

AGJN samples. These studies estimate Gompton-thick 
fractions relative to the overall obscured AGN population 
in the range 30 — 50%, thereby demonstrating that such 
heavily obscured sources represent a sizeable fraction of 
the AGN population in the nearby Universe. Follow- 


in deep X-ray surveys (Norman et al. 

2002 

Tozzi et al. 

2006 Comastri et al. 

2Ulip, an important recent devel- 


opmeht has been the identification of significant samples 


([Brightman & Uedal 

2012 Ceorgantopoulos et al. 

2013 

Buchner et al. 

2U14[ Brightman et al. 2U14|). This has 


ray data, sufficient to constrain the X-ray spectra, along 
with extensive multi-wavelength coverage of X-ray survey 
regions, and new techniques able to determine accurate 


photometric reds hifts for X-ray emitting AGN (Salvato 


et al. 2009 20111. This then offers the exciting possibil¬ 


ity of starting to constrain the evolution of the obscured 
AGN populations, including Compton-thick AGN, pro¬ 
vided that the selection function can be sufficiently well 
understood. 

Further interest in the obscured AGN population lies 
in the nature of the obscuration itself, and its possible 
relationship with other s ource properties. In the stan- 


dard unification picture (Antonucci & Miller 1985 An- 


tonucci 19931, all AGN are surrounded by an optically 
thick toroidal structure relatively close (parsec-scale) to 
the central engine, which can obscure the line of sight de¬ 
pending on the viewing angle. Alternatively, or addition¬ 
ally, obscur ation can occur at galactic scales (Maiolino & 
Rieke|[l9^ I. Observations show that AGN host galaxies 
are massive and lie, at least in a n average sense, on the 
main sequence of star-formation (Santini et al.|2012). At 
moderate redshifts, 2 « 1 — 2, such galaxies are known 
to be gas-rich, wit h gas contents 3 to 10 times larger 
than local samples (Tacconi et al. 20131. It is therefore 
possible that obscuration in moderate redshift AGN is 
associated with the same gas fuelling both star forma- 
tion and the accretion proce ss itself. In some scenarios 


(Hopkins et al. 2006a 20111, obscured AGN represent 


a distinct phase in the co-evolution of the galaxy and 
its central black hole, with energy output from accre¬ 
tion sweeping up gas from the surrounding s with poten¬ 
tially profound effects on star formation (Silk & Rees 
1998 Fabian||1999 King |2003 1. Determining accurately 
the traction of obscured ACJN, including Gompton-thick 
AGN, as a function of other parameters such as lumi¬ 
nosity and redshift is critical for both unification and co¬ 


evolution models. Previous work has provided evidence 
that the obscured fraction depends on luminosity, with 


1 Ueda et al. 

2003 

Treister et al.|2004 

Akylas et al.|2006 

Basinger 21) 


E 

orero et al. 2UU9 Ueda et al. 121) 141), with 


local samples suggesting that the obscured traction may 
peak at around Lx 10"^^ erg/s (Brightman & Nandra 


2011b Burlon et al. 20111. More controversially, it has 


to high redshift (La Franca et al. 20051 

Hasinger 2008 

iirightman k, Uecla|2(J12 

[Iwasawa et al.|2 

OliilIvKo e( al 

20141 but see [Ueda et al. 

20031 [Akylas et al.|2006[[Ebrero 

et al. 12009). 


porting data needed to identify Compton-thick AGN also 
enable an accurate determination of the obscured AGN 
fraction and its evolution, again with the proviso that 
the selection function must be well understood. This 
is, however, very challenging, given the limited photon 
statistics in deep X-ray data, and associated uncertain¬ 
ties in X-ray spectral analysis, combined with additional 
uncertainties e.g. in photometric redshifts. Moreover, 
studies of nearby objects show that the X-ray spectra of 
AGN are complex, including multiple emission compo¬ 
nents, e.g. absorption, reflection or scattering of direct 
AGN emission. 

In this paper we develop a novel non-parametric 
method for determining the space density of AGN as a 
function of accretion luminosity, redshift and hydrogen 
co lumn density. We buil d on the X-ray spectral analysis 


of Buchner et al. (20141, applying their Bayesian spec- 
tral analysis techhique of a realistic, physically motivated 
model to a multi-layered survey determine the luminos¬ 
ity and level of obscuration in a large sample of X-ray 
selected AGN across a wide range of redshifts. Impor¬ 
tant features of our approach are that, first, all sources 
of uncertainty are consistently factored into the analysis 
and secondly, unlike previous studies, no functional form 
is imposed. Instead, we use a smoothness assumption to 
allow the data to determine the space density of AGN 
and its dependence on luminosity, redshift and the level 
of obscuration along the line of sight. The methodology 
is designed to be informative at regions of the parameter 
space where data are sparse and therefore constraints are 
expected to be loose. 

We adopt a cosmology of Hq = 70kms“^Mpc^^ , 
Dm = 0.3 and Uy = 0.7. Solar abundances are as¬ 
sumed. The galactic photo-electric absorption along the 
line of sight direction is modelled with Nh ~ 8.8 x 10^®, 
1.3 X 10®°, 2.7 X 10®° and 2.2 x 10®i cm"® for the 
GDFS, AE GIS-XD, GOSMOS and XMM-XXL fields re¬ 


spectively (Dickey fc Lockman 1990 Stark et al. 1992 
Kalberla et al.|20^). in this work, luminosity [L) always 


refers to the ihtrinsic (absorption-corrected) luminosity 
in the 2 — 10 keV rest-frame energy range. 

2. DATA 

The determination of the obscured fraction of AGN as 
a function of redshift, accretion luminosity requires good 
coverage of the Lx — z plane. We therefore combine X- 
ray survey fields with different characteristics in terms 
of depth and areal coverage. These include th e Chan¬ 
dra Deep Field South (GDFS, |Xue et al. 20111, the All 
Wavelength Extended Groth strip International Survey 
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(AEGIS, |Davis et aL 2007t, the Cosmo logical evolution 
Survey (CUSIVLUS, Scoville et al.||2UU7| and the equato¬ 
rial region of the XJVLM-AAL survey (El: Pierre). 

In our analysis we focus on the region of the CDFS 
which is covered by the 4 Ms Chandra observations and 
the part of the AEGIS field which has been surveyed by 
Chandra to an exposure of 800 ks (AEGIS-XD, Nandra et 
ah, submitted). In the COSMOS field we use the region 
covered by the Chandra observations perform ed between 


to produce event files from the Observation Data Files 
(ODE) using the epchain and EMCHAIN tasks of SAS for 
the EPIC (European Photon Imaging Camera; Striider] 


November 2006 and June 2007 (C-COSMOS, Elvis et ah, 
20091. Table fT] presents information on the individual 
X-ray fields used in this paper. 

2.1. Chandra observations 

The data reduction, source detection and source cata¬ 


logues for t wo of the Chandra fields, the CDFS (Rangel 
et al.||2013|) and the AEGIS-XD (Nandra et ah, subnnt- 
ted). follow the methodology described by Laird et al. 


(20091. Briefly, hot pixels, cosmic ray afterglows and 


times of anomalously high backgrounds are removed to 
produce clean level-2 event files. These are then aligned 
using bright sources and subsequently merged. Images 
and exposure maps are constructed in four energy bands, 
0.5 — 2, 2 — 7, 5 — 7 and 0.5 — 7keV. A candidate 
source list is created using wavdetect at a low signifi¬ 
cance threshold (10“^). Source and background counts 
are then extracted at each candidate source position. 
The source region corresponds to the 70% encircled en¬ 
ergy fraction (EEF) of the point spread function (PSF). 
The background count region is an annulus with inner 
radius 1.5 times the 90% EEF of the PSF and width 
100 pixel of O.Sarcsec in size. For the determination of 
the background, candidate sources are masked out. For 
each candidate source position, the Poisson probability 
that the observed counts are a background fluctuation 
is compute d. Sources are accep ted if that probability is 
< 4 X 10“® (Nandra et al.|2005 1. In this paper we use the 


hard band (2 — 7keV) selected sample. The X-ray sen¬ 
sitivity curves are estimated by extrapolating the back¬ 
ground counts and exposure maps in the 2-7 keV band 
to the limiting flux of a source in the 0.5-1 0 keV energy 
range by ado pting the methods described in |Georgakakis| 
et al.| (|2008|). For the C-COSMOS survey we use the 2- 
7 keV selecte d X-ray source catalogue presented by |Elvis| 
et al. (20091. The number of ^ray sources detected are 
brokeni down by field in Table The corresponding sen¬ 
sitivity curves are shown in Figure 

2.2. XMM observations 

The Chandra surveys above need to be complemented 
by a shallower and wider X-ray field to place constraints 
at the bright end of the luminosity function. For this we 
use the XMM-XXL survey, which consists of 10 ks XMM 
pointings that cover a total area of 50 deg^ split into two 
equal size fields. In this paper we focus on the equato¬ 
rial sub-region of the XMM-XXL. The data reduction, 
source detection and sens itivity map construction follow 
the methods described in Georgakakis & Nandra (20111. 
A full description of the XMM-XXL X-ray source cata- 
logue generation are presented by Liu et al. (in prep.). 
The most salient details of those steps are outlined here. 

The XMM observations were reduced using the Sci¬ 
ence Analysis System (SAS) version 12. The first step is 


et al.|2001 Turner et al. 2001) PN and MOS detectors re- 
spectively. Pixels along the edges of the CCDs of the PN 
and MOS detectors are removed because their inclusion 
often results to spurious detections. Flaring background 
periods are identified and e xcluded using a met hodology 
similar to that described in Nandra et al. (2007). Images 
and exposure maps in celestial coordinates with pixel size 
of 4.35 arcsec are constructed in 5 energy bands, 0.5 — 8, 
0.5 — 2, 2 — 8, 5 —Sand 7.5 —12keV. All overlapping EPIC 
images are merged prior to source detection to increase 
the sensivity to point sources. The detection algorithm 
is applied independently to each of the 5 spectral bands 
defined above. 

The sou rce detection meth odology is similar to that de¬ 
scribed in Laird et al. (20091 in the case of Chandra data. 
Source candidates are identified using the wavelet-based 
EWAVELET source detection task of SAS at a low threshold 
of 4(7 above the background, where a is the RMS of the 
background counts. For each candidate source the Pois¬ 
son probability of a random background fluctuation is 
estimated. This step involves the extraction of the total 
counts at the position of the source and the determina¬ 
tion of the local background value. To match the asym¬ 
metric PSF of XMM, especially off-axis, elliptical aper¬ 
tures were used from the XMM EPIC PSF parametrisa- 
tion of Georgakakis & Nandra (20111. The count extrac¬ 
tion region is obtained by scaling the elliptical apertures 
to contain 70 per cent of the PSF EEF. The total counts 
at a candidate source position, T, is the sum of the ex¬ 
tracted counts from individual EPIC cameras. For each 
source the local background is estimated by first masking 
out all detections within 4 arcmin of the source position 
using an elliptical aperture that corresponds to the 80 per 
cent EEF ellipse. The counts from individual EPIC cam¬ 
eras are then extracted using elliptical annuli centred on 
the source with inner and outer semi-major axes of 5 and 
15 pixels (0.36 and 1.09 arcmin) respectively, while keep¬ 
ing the same shape as the elliptical aperture in terms of 
rotation and ellipticity. The mean local background, B, 
is then estimated by summing up the background counts 
from individual EPIC cameras after scaling them down 
to the area of the source count extraction region. The 
Poisson probability P{T, B) that the extracted counts 
at the source position, T, are a random fluctuation of 
the background is calculated. We consider as sources the 
detections with P(r, R) < 4 x 10“®. The above method¬ 
ology is optimised for the detection of point sources. The 
final catalogue however, includes extended X-ray sources 
associated with hot gas from galaxy clusters or groups. 
Also, the extended X-ray emission regions of bright clus¬ 
ters are often split into multiple spurious detections by 
our source detection pipeline. The emldetect task of 
SAS is used to identify extended sources (i.e. groups or 
clusters) and spurious detections. Point sources for which 
EMLDETECT failed to determine a reliable fit and hence 
are considered spurious and excluded from the analy¬ 
sis. The EPOSCORR task of SAS is used to correct for 
systematic errors in the astrometric positions of X-ray 
sources by cross-correl ating with optically sources in the 
SDSS-DR8 catalogue (Aihara et al. 2011) with magni¬ 
tudes r < 22 mag. 
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TABLE 1 
Survey fields 


Survey 

CDFS 

AEGIS-XD 

C-COSMOS 

XMM-XXL 

Survey area 

464 arcmin^ 

1010 arcmin^ 

0.9 deg^ 

20 deg2 

Total/central exposure time 

4 Ms 

2.4 Ms/800 ks 

1.8 Ms/160 ks 

10 ks 


TABLE 2 

Sample statistics for the individual data extraction steps. 


Section 

Extraction step 

CDFS 

AEGIS-XD 

C-COSMOS 

XMM-XXL 

total 

2.1,2.2 

X-ray hard-band detected 

326 

574 

1016 

206 

2122 

2.3 

Association with optical / IR 

315 

559 

1016 

206 

2096 

2.4 

Redshift information: 

spec-z/photo-z/no-z/removed stars 

180/131/11/4 

227/322/15/10 

491/519/0/6 

174/0/32/0 

1072/986/58/20 

2.5 

X-ray spectral extraction and data analysis: 
successful/failed extraction 

321/1 

564/- 

1010/- 

206/- 

2101/1 

2.6 

galaxies removed (based on 
^2-10 keV < 10'*^ erg/s) 

20 

11 

14 

1 

46 

2.7 

Objects used for LF analysis 

302 

553 

996 

205 

2056 


Total number of sources selected. 


The flux of each source in different spectral bands is 
estimated by assuming a power-law X-ray spectrum with 
r = 1.4, i.e. similar to the XRB, absorbed by the appro¬ 
priate Galactic hydrogen co lumn density. The la tter is 
derived from the HI map of Kalberla et al. (20051 using 
the right ascension and declination ot the aimpoint of 
each XMM observation and the NH task of ftools. The 
energy to flux conversion factors are such that the counts 
from the 0.5-2, 0.5-8, 2-8, 5-8 and 7.5-12 keV bands are 
transformed to fluxes in the 0.5-2, 0.5-10, 2-10, 5-10 and 
7.5-12 keV bands respectively. 

In this work we use the XMM-XXL point source 
sub-sample that is selected in the 2-8 keV band. 
We minimise optical identification and spectro¬ 
scopic redshift determination incompleteness (see 
next sections) by applying a bright flux cut, 
/x(2 — lOkeV) > 7 X 10“^^ ergs“^ cm“^. The total 
number of XMM-XXL X-ray sources used in the 
analysis is 206 (see Table |^. 

The construc tion of the sensitivity ma ps follows the 
methodology of Georgakakis et al. (20081. The 2-10 keV 
sensitivity curve ot the sample is plotted in Figure 
The window function of the spectroscopic follow- 
up observations in the XMM-XXL held (see Section 
\2A\ is taken into account in this calculation. Also, 
in addition to the Poisson false detection probabil¬ 
ity threshold of < 4 x 10“® in the sensitivity map 
calculation we also take into account the flux cut 
/x(2 — lOkeV) > 7 X 10“^'^ergs“^ cm“^, by estimating 
at each survey position the probability of measuring a 
flux above this limit. 

The number of X-ray sources detected - about 2000 
in total - are broken down by held in Table The 
corresponding sensitivity curves are shown in Figure 
We associate these detected X-ray source positions to 
optical counterparts in order to obtain redshifts. 

2.3. Association with optical/IR counterparts 

The identification of the X-ray sources with optical 
or infrared counterparts in the AEGIS-XD, GOSMOS 
and XMM-XXL used the Likelihood Ratio method of 


association of X-ray sources with optical/infrared coun¬ 
terparts are presented by Nandra et al. (snbmitted). 
They used the multi-waveband photometric catalogue 
provided by the Rainbow Gosmological Surveys Pa t abase 
(Perez-Gonzalez et al. 2008 Barro et al. ]|2011a|b| . 


'The cou nterparts o: 


G- 0(JSlV l(Jb X-ray sources are 


taken from Givano et al. (|2012|). For the identi fication 
they 


et al. 


used the 7-band selected optical sample of 


et al. 


et al. 


Gapakl 


,pa 

20071, and the AT-band photometry of|McGracken| 


and the IRAG-3.6^m catalogue of banders 


X-ray sources in the XMM-XXL survey were matched 
to the SDSS-DR8 photometric catalogue [Alhara^e^aL 


(|2011t following the methods described in Georgakakis 
fcblaiidr^ ([20 ^. 


IHsu et al. (|2bl4l presents the counterparts of the 
GDFS X-ray sources. They used a sophisticated 
Bayesian version of the Likelihood Ratio m ethod which 
is based on t he probabilistic formalism of [Budavari fc| 
Szalay (20081. The photometric catalogues used to iden¬ 
tity xAay sources include the G ANDELS iJ-band se¬ 
lected multi-wavelength catalog in Guo et al. J 201 3 1, the 
MUSYG catalogue presented by Gardamone et al. 12010 1 
and the TE NIS near-infrared s Gected source catalogue 
described by Hsieh et al. (2012|. 


2.4. Redshift estimation 

The X-ray survey fields used in this paper benefit 
from extensive spectroscopic campaigns that also tar¬ 
get specifically X-ray sources. In the GDFS, we used 
the spectroscopic re dshifts compiled by N. Hathi (private 
communication, see|Hsu et al.||2014|). Spectroscopic red¬ 
shift measurements ot X-ray sources in the AEGIS held 
are extracted from the compilation presented in Nan- 
dra et al. (submitted ) which included also the DEEP2 
(Newman et al.||2012t and DEEP3 galaxy redshift sur- 


veys (|Gooper et al.||201l| 
carried out at the iVliVl'i' i 


20121 as well as observations 
using the Hectospec fibre spec- 


Sutherland & Saunders (19921. Specific details on the 


trograph (Goil et al.|2009). Redshifts in the G-GOSMOS 
are used irom the compilation of Givano et al. (20121 
which includes the pu blic releases ol th e VllVlUS/zGUS- 
MOS bright project (Lilly et al. 20091 and the Magel- 
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In the case of XMM-X 

.XL, optical spectroscopy is 
the Baryon Oscillation Spec- 

from Stalin et al. 

(20101, 

troscopic Survey i 

roSS; 

Dawson et al. 2013| 

Bolton 

et al. 2012 Smee et al. 

|2013p as well as dec 

heated 


tSloan Digit al tSky Purvey 


Gunn et al. 2006 SDSS-III: Eisenstein et al. 2011 
ancillary science observations, which targeted specib- 
cally X-ray sources in the equatorial XMM-XXL field 
(PI: A. Merloni, A. Georgakakis). Targets were se¬ 
lected to have /x(0.5 — lOkeV) > 10“^^ergs“^ cin“^ 
and 17 < r < 22.5, where r corresponds either 
the PSF magnitude in the case of optical unresolved 
sources (SDSS type=6) or the model magnitude for re¬ 
solved sources (see Menzel et al., in prep.). At the 
flux limit fxi‘2 — lOkeV) > 7 x 10“^'^ ergs“^ cm“^, 84% 
(174/207) of the XMM-XXL sources have secure redshift 
measurements. 

The GDPS, AEGIS-XD and GOSMOS fields have mul¬ 
tiwavelength photometric observations that allow pho¬ 
tometric redshift estimates for sou rces that lack spec¬ 


trosco py. The methods developed by Salvato et al. (2009 


20111 are adopted to achieve photometric redshift ac- 


curacies for X-ray AGN comparable to galaxy samples. 
The photometric redshifts catalogues are presented in 


Hsu et al.| (2014) for the GDPS, N andra et al. (submit¬ 
ted) tor tne AEGIS-XD held and (Salvato et al.|2011l 


for the GOSMOS held. The photometric redshiits are 
included in the analysis in the form of probability distri¬ 
bution functions, incorporating systematic uncertainties 
(see below). A by-product of the photometric redshift 
determination is the characterisation of the Spectral En¬ 
ergy Distribution (SED) of X-ray sources. This informa¬ 
tion is used to identify candidate Galactic stars in the 
X-ray sample and ex clude from the analysis. Following 


Salvato et al. (20111 if stellar templates provide an im- 
proved ht to tne SED of a source, as measured by the 
reduced Xstar < Xgai/l-^i and the source is point¬ 

like in the optical images, then it is considered to be a 
Galactic star candidate. The number of removed stars is 
indicated in Table |2] 

For the XMM-XXL held, multiwavelength photometry 
from the UV to the infrared which homogeneously covers 
the surveyed area is not available. This is essential for 
reliable photometric redshifts, especially in the case of 


bright X-ra y samples like the XMM-XXL (e.g. Salvato 
et al.|[MTT). Therefore, for the number of X-ray sources 


without spectroscopic redshift measurement in that held 
(32/206) we choose not to determine photometric red¬ 
shifts. These sources are still included in the analysis 
by assigning them a flat redshift prior (see below for de¬ 
tails). In the XMM-XXL bright sub-sample used here, 
there are no stars based on inspection of the images. 

To summarise, the redshift determination falls in one 
of three cases for each X-ray source: 


1. No redshift information available (58 sources): 
This is the case if no spectroscopic redshifts are 
available, and photometric redshift estimation is 
not possible due to limited photometry. This case 
also applies when no secure association was found. 
Such a source is associated with a flat redshift dis¬ 
tribution in the interval z = 0.001 — 7. 


2. Spectroscopic redshifts (1072 sources): Wherever 
secure spectroscopic redshifts are available, they 
are used directly, i.e. do not associate uncertainty 
to them. Different surveys use different conven¬ 
tions to define the reliability of the redshift mea¬ 
surement. In this paper we only consider spectro¬ 
scopic redshifts from the top two quality classes of 
any study, which typically corresponds to a proba¬ 
bility better than 95 per cent of being correct. 

3. Photometric redshifts (986 sources): These are in¬ 
cluded in the analysis in the form of probability 
distribution functions. We also incorporate system¬ 
atic uncertainties e.g. due to incorrect association 
to an optical counterpart (see Appendix B). 

Table shows a breakdown of how many sources fall into 
each category. The distribution of redshifts is plotted in 
Figure 


2.5. Extraction and analysis of X-ray spectra 

For the Ghan dra surveys, t he ACIS EXTRACT (AE) soft¬ 
ware package (Broos et al. 20101 was used to extract 
spectra for each source. 'I' he extraction followed the same 
methodology described in |Brightman et al. 

each p 


tially, each source and each pointing is dea. 


(20141. Ini- 
t with sepa¬ 


rately. AE simulates the PSFs at each source position. 
Regions enclosing 90% PSF at 1.5keV were used to ex¬ 
tract source spectra. The background regions are con¬ 
structed around the sources such that they contain at 
least 100 counts, with other sources masked out. AE 
also constructs local response matrix flies (RMF) and 
auxiliary matrix flies (ARF) using ray-tracing. As a fi¬ 
nal step, AE merges the extracted spectra so that each 
source has a single source spectrum, a single local back¬ 
ground spectrum, ARF and RMF. 

In the XMM-XXL held. X-ray source and background 
spectra were extracted separately for each EPIG cam¬ 
era using SAS version 12.0.1. The source extraction re¬ 
gions were chosen to maximise the signal-to-noise ratio 
using the eregionanalyse task of SAS. In the case of 
EPIG-PN camera, the background region has a radius 
of 90 arcsec and is placed, after visual inspection, on a 
source-free region on the same GGD and raw-Y axis as 
the corresponding source. In the case of EPIG-MOS be¬ 
cause of the larger size of individual GGDs it is possible 
to define the background region as an annulus with the 
inner and outer radius were 40, 110 arcsec respectively. 
If that was not possible, because e.g. a source might lie 
close to GGD gaps then the background regions was a 
circle of 90 arcsec radius placed, on a source-free region 
on the same GGD as the corresponding source. ARFs 
and RMFs are generated using the ARFGEN and RMFGEN 
tasks of SAS. Finally, the ARFs, RMFs, source and back¬ 
ground spectra of the same object observed by different 
EPIG cameras are merged by weighting with the expo¬ 
sure time to produce unique spectral products for each 
source. 

Each X-ray spectrum was analysed with spectral 
models fo ll owing the Bayesian methodology of |Buch-| 
ner et al. (20141. We only consider their best model, 


torus+pexmon+scattering, which consists of 

1. an intrinsic power law spectrum modified by photo¬ 
electric absorption and Gompton scattering from 
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Fig. 1.— Luminosity-redshift plot of the full sample. The intrinsic 2 — 10 keV luminosity in erg/s and redshift is plotted. For this 
visualisation, we colour sources based on logNfj (in units of cm~^ as follows: 20 — 22, “Unobscured”, red; 22 — 24, “Compton-thin Obscured”, 
blue; 24 — 26, “Compton-thick”, green). If more than 50% of the Nfj posterior probability lies within one of the intervals above the source 
is colour-coded accordingly. Due to the heavy suppression of the flux of Compton thick objects by absorption, they are typically only 
detectable at higher intrinsic luminosity compared to unobscured or Compton-thin sources. Some objects lie below the L = erg/s 
limit and are not used in this work. The inset shows the Nu histogram of our sample, while the top and left axes show histograms of the 
redshift and luminosity distributions in gray. These plots are constructed by drawing a random posterior sample for each object, and 
thus includes the uncertainty in the parameter estimation as well as the intrinsic distribution. Sources above redshift z = Z which have 
spectroscopic redshift estimates are indicated with large black circles. 


Sensitivity 


10^ 

10 ° 


1 

S’ 10'^ 

■a 

03 

> 10 -^ 

< 

10 '^ 


10-17 ^q-16 ^q-15 ^q-14 ^q-13 ^q-12 

Flux [erg/s] (2-10 keV) 

Fig. 2. — Area curve. The bottom plot shows all X-ray sensitiv¬ 
ity curves on a logarithmic scale for the sum of all fields (dotted 
magenta line) and the individual fields. The top plot shows the 
sensitivity curves for the individual fields on a linear scale and 
normalised to their respective maximum area (see Table to in¬ 
dicate the flux limit for detection. The XMM-XXL curve is limited 
in flux by 7 X 10“^^ erg/s, but we incorporate the uncertainty of 
measuring a higher flux due to Poisson variance. This introduces 
a smooth transition. 

an obscurer with toroidal geometry simulated by 



the torus model of Brightman & Nandra (2011al, 

2. a Compton-reflecti on component appro ximated by 
the pexmon model (Nandra et al.|[2ro7l and 


3. a soft scattering component, which is parame- 
terised by an unabsorbed power law with the same 
spectral slope as the intrinsic one and normalisa¬ 
tion that cannot exceed 10% of the intrinsic power- 
law spectrum. 


This model has five free parameters, the slope of the in¬ 
trinsic power-law spectrum, T, the line-of-sight column 
density of the obscurer, Nh, the normalisations of the 
intrinsic power-law, Compton-reflection an d soft scatter¬ 
ing cornponents . For the torus model of [Brightman fc 


Nandra (2011al we fix the opening angle parameter to 
45” and the viewing angle parameter to the maximum 
allowed, 87°, i.e. edge on. For pexmon component we 
use the same incident power-law photon index as the in¬ 
trinsic one spectrum and fix the inclination to 60°. The 
normalisation of the pexmon (i? parameter) is modelled 
relative to that of the torus component and is allowed 
to vary between 0.1 and 100. In the analysis, we use flat 
priors on the luminosity/normalisation parameters and 
the column density, and an informed Gaussian prior of 
mean 1. 95 and standard deviatio n 0.15 for the photon 
index F (Nandra & Pounds||1994 ). 

In the spectral analysis, redshifts are either fixed, when 
spectroscopy is available, or are included in the form 
of priors that follow a probability distribution function 
(photo metr ic redshifts and missing counter parts, see 
Section 2.4 above). 

Following the methodology of Buchner et al. (20141, we 
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use the 0.5—8 keV spectrum and the redshift information 
to infer posterior distributions for the X-ray luminosity 
Lx, redshift z and column density Nh- 

2.6. Sample selection 

Based on the initial hard-band detection, normal 
galaxies (i.e. non-AGN) should already be largely ab¬ 
sent from our sample, except for strongly star-bursting 
galaxies. To exclude these in a conservative manner, 
we remove intrinsically faint sources that have Lx < 
10'^^ erg/s with 90% probability. The number of objects 
used in the luminosity function analysis is shown in Ta¬ 
ble Figure [J shows how the sample is spread in the 
luminosity-reddrift plane. 

3. METHODOLOGY 

In this work, we want to study the evolution of AGN 
obscuration. More specifically, we are interested in the 
distribution of the population in X-ray luminosity L and 
column density Nh, i.e. the luminosity function. Our 
sample is a specific draw from this distribution, with a 
sampling bias in detection and incomplete information on 
each object. Inference meth ods for such scenario s have 
been known for a long time (Marshall et al.|[r983 1 . 


3.1. Luminosity function analysis 

Our general approach to establish the evolution of the 
obscured and unobscured populations of AGN is to de¬ 
termine the X-ray luminosity function, accounting for 
all measurement uncertainties, which can then be anal¬ 
ysed as a function of obscuration to determine the abso¬ 
lute and relative number densities of obscured and un¬ 
obscured sources. 

Direct visualisation of the data is difficult if we want to 
stay true to the uncertainties. Faint objects, which dom¬ 
inate the sample, are highly uncertain in their proper¬ 
ties, namely the intrinsic luminosity L, obscuration Nh, 
but also their redshift z in the case of photometric es¬ 
timates. For instance, Gompton-thick AGN can have 
large “probability clouds” for their parameters. This pro¬ 
hibits us from assigning objects to bins for visualisation. 
For direct visualisation, three approaches can be consid¬ 
ered: (1) assigning each object to a random luminosity 
bin based on its probability distribution, and then esti¬ 
mating the density in the bin, (2) assigning each object 
to each luminosity bin with a probability weight, and 
then estimating the density in the bin, (3) computing 
for each luminosity bin the number of objects that have 
a higher luminosity with e.g. 90% probability. Method 
(2) has the difficulty that the “number” in each bin is 
no longer integer - requiring interpolation of the Poisson 
distribution formula. The methods (1) and (2) assume 
a frequency interpretation of the uncertainty probability 
distributions - which is not reasonable as every object is 
different and the sample size is small. Method (3) may 
be useful for checking whether data and model agree, 
but does not yield an intuitive visualisation. We thus 
find none of these methods satisfying, and develop a new 
approach (next section), which visualises the data and 
estimates the luminosity function at the same time. 

In the Appendix A, we review the statistical footing 
of analysing population demographics by re viewing and 
combi ning the works of Loredo (2004) and Kelly et al. 
(20081. In this work, we use the usual Poisson likelihood: 


UkI^-Pidk,D\C)-^dC 

-^(C) • 


( 1 ) 


where C = {logL, z, log Nh}, p{dk, D\C) represents the 
results of the spectral analysis of data dk from the de¬ 
tected object k, which is weighted by the luminosity 
function (j). The integral in the denominator computes 
the expected number of sources by convolving the lu¬ 
minosity function with the sensitivity curve A (shown 
in Figure 1^ as a dotted magenta line). The conver¬ 
sion from 77, z, Nh to sensitive area requires a spectral 
model, for which we use the torus model, and we aver¬ 
age over F in exact correspondence to the Gamma prior 
used in our analysis. Here, the additional sca tter ing com¬ 
ponents (+pexmon+scattering, see Section 2.5) are not 
used. Both contribute on ly minimally to the 2 — 10 keV 
flux: [Buchner et al.|20T4 presents an extreme example in 
detail, where both components show the strongest con¬ 
tribution in their sample, but the 2—10 keV flux is only 
affected by 10% (0.04 dex). Future studies may take 
these components and their normalisations into account 
as well. The Appendix A gives an extensive derivation 
and caveats for equation ^ 


3.2. Non-parametric approach 

We would like to use the power and safety of a 
likelihood-based analysis but without the rigidity of a 
functional form, allowing discovery of the shape of the 
luminosity function. Our method for analysing the lu¬ 
minosity function is thus, in a simplistic description, to 
fit a three-dimensional (L, z, Nh) histogram as the lu¬ 
minosity function model. 

Using e.g. 10 x 10 x 10 = 1000 bins already means 
that the problem is largely under-determined. We thus 
need to input additional knowledge. Here we make the 
reasonable assumption that the function does not vary 
rapidly between neighbouring bins (in particular for the 
redshift). We encode this smoothness prior in two ap¬ 
proaches, both using the Normal distribution to penalise 
large deviations. These two priors, “constant-value” and 
“constant-slope” are explained below and illustrated in 
Figure 

• “constant-value”: This prior retains the current 

value unless constraints are imposed by the data. 
The value of a bin should scatter around its neigh¬ 
bour density value = Normal(log ())i, cr) 

with an allowed correlation width a for each axis 
(ctl, cr^). 

• “constant-slope”: This prior keeps power law slopes 
intact unless constraints are imposed by the data. 
The log-density slope between bins should scatter 
around its neighbours slope log 0^+2 = log(/>i_|_i -I- 
Normal(log— log^i, cr), with the deviation 
from the slope for each axis {cjl, ctz)- 

For the Nh dimension, we always use the constant- 
value prior, as in our view a power law dependence is not 
appropriate here. With this simple prescription, we can 
recover a smooth field by fitting a model, whose shape is 
driven by the data. 
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prior: constant-value 


prior: constant-slope 



1 2 3 


2 3 


bln 


bin 


Fig. 3. — Illustration of the two smoothness priors used. The 
left panel shows the “constant-value” prior: The extrapolation from 
well-constrained step-function bins (black) to a neighbouring bin 
(red) is done by assuming the same value with a fixed uncertainty a, 
whose value encodes the assumed correlation strength or smooth¬ 
ness. The right panel illustrates the “constant-slope” prior. The 
value for the neighbouring bin (red) is predicted by continuing the 
same slope from the black points. As the prior is defined in logarith¬ 
mic units of the density, this behaviour corresponds to preferring 
a powerlaw. 


nique to ensure rapid mixing of the Markov chain by 
avoiding cyclic explorations. Our Stan model is shown 
in Algorithm^ in the Appendix. 

4. RESULTS 

A few words on the visualisation are warranted. In 
all relevant figures we plot the posterior distributions of 
our three-dimensional step function model from various 
axis views {L, z and Nh)- We always show the median 
result for the constant-value prior as a dashed line, and 
the median result for the constant-slope prior as a solid 
line. A feature of our methodology is that uncertainties 
are realistic and reflect regions of parameter space where 
data are sparse. Whenever the data constrain the result 
well, the results from either prior prescription (constant- 
value or constant-slope) will be the same. Only when 
constraints from the data are poor, the results differ. 
We thus take the difference in the reconstructions as a 
indication of whether the data or the priors dominate the 
result. We plot the 10%-90% quantile hatched regions 
as a measure of the uncertainty, by taking together the 
posterior samples from both priors. 


Whenever the data constrain the result well, the re¬ 
sults from either prior prescription (constant-value and 
constant-slope) will be the same. Where the constraints 
are poor, the results will differ depending on the adopted 
prior. We thus take the difference in the reconstructions 
as an indication of whether the data or the priors domi¬ 
nate, and in the latter case the difference between the two 
is an indication of the uncertainty in the determination. 

It should be stressed that the choice of binning and 
in particular the correlation strengths ap, <7 Nh can 
influence the result. Motivated by the number of data 
points in each bin, we choose our pixelation as 11 bins of 
logarithmically spaced luminosity 42 — 46 (units erg/s), 
redshift bin edges 0.001, 0.1, 0.3, 0.5, 0.75, 1, 1.25, 1.5, 
1.8, 2.1, 2.7, 3.2, 4, 7 and log Nh bin edges 20, 21, 22, 
23, 23.5, 24, 26. The correlation strengths a is defined 
for neighbouring bins, and their choice is important. No¬ 
tably, if cr is too small, the model will be flattened out, as 
the prior dominates, while if a is too large, no smooth¬ 
ness assumption is used, and uncertainties will be large. 
Ideally, we would like to recover a from the data them¬ 
selves. Unfortunately, our tests show that this computa¬ 
tion is not numerically stable. We do And however that 
above some value of a, the results are stable regardless 
of the choice of cr. We thus just choose reasonable values 
for cr, namely up = 0.5, cr^ = 0.5 and (JNh — 0.75. This 
encodes, roughly speaking, that neighbouring bins have 
the same order of magnitude in space density. These val¬ 
ues have been chosen after a few initial tests, but were 
not tuned to give optimal results. Rather, we believe, 
they are one possible, and reasonable, a priori choice. 

We are also highly interested in the uncertainties of our 
smooth field reconstruction method. Bins with data will 
be tightly constrained, while bins without information 
will have increasing uncertainty with distance from the 
data. Making uncertainty estimates with so many pa¬ 
rameters is not trivial, especially as the parameters are 
correlated by definition. We use a Hamil tonian Markov 
Chain Monte Carlo code named “Stan” (Stan Develop- 


ment Team| 20141. St an uses the sophisticated INo-G^ 


4.1. Total luminosity function 

Figure shows our non-parametric total X-ray lumi¬ 
nosity function, integrated along the Nh axis. Large 
uncertainties are visible at the bright end at low redshift 
(z < 0.3), where our sample is small due to the lim¬ 
ited cosmological volume. At high redshift (z > 2), the 
faint end becomes uncertain as low X-ray fluxes limit the 
number of detected sources. 

Our method captures the general shape, normalisation 
and evolution with redshift of the X-ray luminosity func¬ 
tion (XLF) as inferred or assumed in previous parametric 
determinations. In Figure we show in particular the 
comparis on with a r e cent comprehensive study of the 
XLF by Ueda et al. (2014|). The overall shape of the 
luminosity function is a double power-law with a break 
or bend at a characteristic luminosity (L*), the value 
of which increases with redshift. As found in previous 
studies, the space density shows a rapid evolution up to 
around z ^ 1 at all luminosities, being most prominent 
at high luminosities due to the positive evolution of L*. 
We fi nd that the posit ive evolution continues up to z ~ 3 


(as in Aird et al.|2010 1 above which the population starts 


to decline. It is important to emphasise that in our non- 
parametric approach these features are imposed by the 
data and not by any assumptions about the functional 
form of the X-ray luminosity function. 

While the general behaviour of the population is simi¬ 
lar, Figure^ also shows important differences compared 
to some previous parametric studies, specifically that of 


Ueda et al.|(|2014 1. At the highest redshift bin (z = 4—7), 
the space density drops sharply tow ards high redshifts 


in their XLF reconstruction (see also Ka lfountzou et al.| 
2014 Vito et al.|2014 Civano et al.|201l| . However, our 


Turn Sampler (NUTS, Hoffman & Gelman 20111 tech- 


reconstruction remains, as suggested by the priors, at 
space densities comparable to the previous redshift bin 
z = 3 — 4. In our data, there seems to be no evidence of 
a steep decline with redshift. The difference may be due 
to the large uncertainties in the redshift est imates used 
i n this work. At intermediate redshifts, the |Ueda et al.| 
(20141 XLF shows a sharp flattening below a luminosity 
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1^^ This work: 90% credible interval 

This work: using constant-value prior 
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Fig. 4. — Total X-ray luminosity function in the 2 — 10 keV spectral band. Each panel corresponds to different redshift intervals. We 
show our field reconstruction based on a step function in black. The dashed black lines refer to the prior preferring to keep the space 
density constant, while the solid black line prefers to keep the slope of the density (in z and L) constant. Whenever the data constrain 
the result well, the space density estimate from either prior prescription (constant-value or constant-slope) is the same. Only when they 
extrapolate away from any constraints, the results differ. We thus take the difference in the reconstructions as an indication of whether 
the data or the priors dominate the result. The hatched regions indicate our measure of the uncertai nty, using the 10%- 9Q% quantiles of 
the posterior samples from both priors together. The orange thin solid line shows the reconstruction by|Ueda et al.|||2014|l. The dotted red 
curve is their local {z = 0.1) luminosity funct ion kept constant across all panels for comparison. There are important differences between 
our luminosity function estimates and that of |Ueda et al. pjni discussed in the text. 


TABLE 3 

Key statistics on the fraction of obscured and Compton-thick AGN. 



Cosmic time average^ 

At 2 = 0®’'" 

Maximum^’^ 

C 

'^Maximum 

Obscured fraction (> 10^^ cm ^): 

75fl% 

77l^% 

83l®% 

> 2.25 

Compton-thick fraction (> 10^^ cm~^): 

38t*% 

39tl% 

46l®% 

(unconstrained) 


^ These fractions relative to the total space density of the population are estimated by integrating the X-ray luminosity 
function over cosmic time and within X-ray luminosity range L(2 — 10 keV) = — lO"^® erg/s. The uncertainties 

are computed using the 10% and 90% quantiles of the posterior distribution, i.e. the true value is bracketed with 90% 
probability. 

^ For this estimate we used our lowest redshift bin. 

^ The peak is computed by identifying the maximum value (fraction) and location (redshift) in each posterior realisation, 
and considering the distribution of each series. This leads only an upper limit for the peak location for obscured AGN. 
For Compton-thick AGN, no upper or lower limit could be determined. 


of around 10^^ erg/s, after which it steepens again at the 
lowest luminosities probed by the study. This behaviour 
is most apparent at redshifts z ~ 1 — 2 and is only easy 
to parameterise using the Luminos i ty De p endent Density 


Evolution model (e.g. Ueda et al. 

(20031; Hasinger et al. 

(20051; 

yUverman et al.|(2UU8p). 

Uur analysis does not 

support such a behaviour. This is most apparent in the 


redshift range z = 0.8 — 2.1, where our non-parametric 
field analysis requires a significantly larger space density 
of AGN in the critical luminosity range of erg/s. 


Furthermore we see no strong evidence of a change in 
the faint-end slope in the individual panels of Figure 
This may bring into question the need for LDDE to ex¬ 
plain the form and evolution of the XLF. As our focus in 
this paper is on the demographics and evolution of ob¬ 
scuration, rather the luminosity function per se, we defer 
further discussion of this point to later work. 


4.2. Obscured and Compton-thick fractions 
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We now report on the fraction of obscured AGN by 
comparing the space density above Nh > 10^^ cm“^ to 
the total. Our non-parametric approach allows us to ex¬ 
plore this fraction in a model independent way by inte¬ 
grating the luminosity function over luminosity and cos¬ 
mic time. We computed these fractions using only the 
luminosity range L = 10^3.2-46 erg/s, i.e. without the 
lowest three luminosity bins. The choice of the luminos¬ 
ity range for the presentation of the results is to minimise 
uncertainties associated with the typically looser con¬ 
straints our methodology yields at the faint-end of the 
X-ray luminosity function. At virtually all redshifts the 
AGN space density is better determined at luminosities 
L > erg/s. In Table we find that the fraction 

of obscured objects in the Universe is 75%, with nar¬ 
row uncertainties. The fraction of Gompton-thick AGN 
{Nh > cm“^) is approximately 35%. 

Additionally, we test whether these fractions are con¬ 
stant over cosmic time. This is done by noting the peak 
fraction and its redshift in each posterior realisation of 
our reconstruction, and averaging the results. As indi¬ 
cated in Table the obscured fraction tends to a higher 
value of almost 85% at z > 2. This shows that the frac¬ 
tion of obscured AGN varies through cosmic time. In 
the following section, we investigate this evolution of the 
obscured fraction further. For the Gompton-thick frac¬ 
tion, no peak can be identified as the redshift is uncon¬ 
strained. This indicates that the Gompton-thick fraction 
is constant over cosmic time at approximately 35%. 

4.3. Obscuration-dependent evolution 

To explore the evolution of the obscured AGN frac¬ 
tion further. Figure plots the evolution of the space 
density of unobscurea {Nh = cm“^), moder¬ 
ately obscured {Nh = cm“^) and Gompton- 

thick {Nh = cm“^) AGN with luminosities 

L > 1043-2 erg/s. In the left panel, we find that mod¬ 
erately obscured and unobscured AGN follow similar 
evolutionary patterns, namely an increase from z = 0 
to z = 1.2 where their space density peaks and a de¬ 
cline at higher redshifts. However, there are also dif¬ 
ferences. Moderately obscured AGN evolve much faster 
from z ~ 0.5 to z ~ 1.5. 

The right panel of Figure shows the space density 
of Gompton-thick AGN. The evolution of the Gompton- 
thick AGN population has larger uncertainties than un¬ 
obscured and moderately obscured Gompton-thin AGN. 
Nevertheless we find that their space density shows a 
broad plateau at z « 1 — 3, followed by a decline to both 
lower and higher redshifts. The evolution of Gompton- 
thick AGN also appears weaker than that of moderately 
obscured Gompton-thin ones. This behaviour is con¬ 
trary to the Nh smoothness prior, which prefers neigh¬ 
bouring Nh bins to have the same value. Thus it can 
be concluded that the data drive the result of different 
evolutions for different obscurations: A strong obscured 
Gompton-thin evolution compared to a weaker evolution 
in both Gompton-thick and unobscured AGN. 

The different evolution of AGN with different levels of 
obscuration is further demonstrated in Figure 1^ where 
the AGN sample is split into finer Nh bins. ^11 sub¬ 
populations experience the same evolutionary pattern, a 
rise from redshift z « 0.5, a peak at z « 1.5 and a decline 


at higher redshift. The AGN that undergo the strongest 
evolution are those with columns densities around Nh = 
1022“23.5 Both unobscured and Gompton-thick 

AGN evolve less strongly. 

4.4. Luminosity-dependence and evolution of the 
obscured fraction 

Previous studies suggest that the fraction o f obscured 
AGN is a function of accretion lumino shy (e.g. Lawrence 


19911|Ueda et al.|2003[|Simpson|20051|Akyias et al.|2UUt)l 

■ ^|2oo8'['~ 


Silverman et al. | |2008[^urlon et al.J|2011[ |Ueda et al.| 

2014||. The luminosity-dependence of the obscured frac¬ 
tion is thought to be related to the nature of the ob¬ 
scurer, e.g. reducing it in physical extent. In this sec¬ 
tion we investigate this issue using our model indepen¬ 
dent approach. To illuminate the behaviour of the ob¬ 
scured population we analyse the behaviour of Gompton- 
thin and Gompton-thick AGN separately. This section 
only considers the obscured AGN that are Gompton-thin 
{Nh < 1024 cm“2). To this end, we define a new quan¬ 
tity, namely the obscured fraction of the Gompton-thin 
AGN (Gompton-thin obscured fraction, GTNOF) 


TNOF ■= = 10"-" cm-2] 

■ (I)[Nh = 1020-24 cm-2] • 


( 2 ) 


In the top row of Figure we find that the GTNOF is 
a strong function of luminosity. There is evidence for a 
peak at a certain luminosity and decline at both brighter 
and fainter luminosities (the red star symbol provides a 
constant reference point). Interestingly, the luminosity 
where the obscured AGN fraction peaks appears to be a 
function of redshift. With increasing redshift, the drop 
of the GTNOF at bright luminosities occurs at higher 
luminosities. These r esults are in some agreement with 
the recent analysis in Ueda et al. (20141. One difference 
however is at the brightest luminosities, where our non- 
parametric method tends towards a higher value (5 0%). 
Similar high obsc ured fractions were sug gested by |Iwa-| 


sawa et al. 


(20121 and Vito et al. (| 2014L as opp osed to 


20% in tEe local Universe Burlon et al. (2011). Ueda 


et al. (2014) have better statistics at bright iuminosi- 


ties compared to our work because they included more 
wide-area and shallow survey fields in the analysis. It is 
therefore likely that our sample has poor statistics at the 
brightest luminosities which makes the preference of the 
Nh prior towards equipartition apparent. The strength 
of our data is rather at the faint-end of the luminos¬ 
ity function for moderate and high redshifts. There we 
find a significant turnover at low luminosities (right top 
panel of Figure [^. This result is independent of the 
adopted form of the prior, indicating that this behaviour 
is strongly imposed by the data. A similar behaviour 
of a peak lumi nosity and a turnover have been found in 
local s amples (Burlon et al. 2011J_ Brightman & Nandra| 
'2011b[| (top left panel of h'igure corrected assuming a 


constant 20% Gompton-thick fraction in their sample). 
The position of the peak they find is consistent with our 
results at low redshift (red star symbol, L = lO^® erg/s). 

When considering the obscured fraction at L = 
1044 erg/s, the evolution of the peak is imprinted as a 
rise with increasing redshift. This is shown in Figure 
with the same data points from literature. Observers 
considering mostly AGN with luminosity L > lO^^ erg/s 
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Fig. 5.— Redshift evolution of the space density of AGN with Lx > erg/s, for various column densities. Our step function 

reconstruction is represented by points at the bin centre, which are connected by lines (dashed for the constant-value prior, solid lines for 
the constant-slope prior). In the left panel, distinct evolutions for the Compton-thin obscured (blue shaded region, top) and unobscured 
(red shaded region, bottom) can be observed. The right panel plots the evolution of Compton-thick AGN as a green shaded region. To 
facilitate the comparison we also plot the evolution of the unobscured and obscured Compton-thin AGN reconstruction in the case of the 
constant-slope prior (solid lines). All AGN sub-populations split by the level of obscuration experience similar space density evolution, 
which can be described by a rise from 2 = 0.5 to z = 1.25, a broad plateau at 2 = 1.25 — 2.1 and a decline at higher redshift. There is also 
evidence that moderately obscured, Compton-thin AGN {N^j = 10^^ — 10^“^ cm~^) are evolving faster in the redshift interval 0.5 — 4 in 
the sense that they reach peak space densities higher than the other AGN sub-populations. The space density of Compton-thick AGN has 
the highest uncertainty, due to poor statistics in the low-luminosity range (L < 10'^^ erg/s). Nevertheless, there is tentative evidence that 
the evolution of Compton-thick AGN is weaker than that of the Compton-thin obscured AGN (blue), and in fact closer to the evolution of 
the unobscured AGN (bottom red solid line). 
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Fig. 6.— Redshift evolution of space density of AGN split by the level of obscuration. Different panels correspond to different hydrogen 
column density interval as indicated at the top. The plot highlights that the redshift evolution of the AGN space density is strongest for 
the column density bins Nh = 10^^ — 10^® cm~^ and Nh = 10® — 10^^ ® cm“^. 


will see an increase in the fraction of obscured AGN with 
redshift as shown here, while observers considering in¬ 
trinsically faint AGN (L < 10"^^ erg/s) would observe 
the opposite trend (see Figure]^. 

4.5. Evolution of Compton-thick AGN 

Gompton-thick objects have been hypothesised to play 
a major role in the accretion phase of AGN. The iden¬ 
tification of such sources in current X-ray surveys by 
XMM-Newton and Ghandra is challenging and there¬ 
fore their space density as a function of redshift has re¬ 
mained controversial. It is therefore important to place 
constraints on their space density of these sources in the 
context of previous studies. Figurej^shows the evolution 


of the space density evolution of Gompton-thick sources 
inferred by our non-parametric methodology. Three lu¬ 
minosity cuts are show n which allow direct comparison 


with previous st udies (Fiore et al. 2008 2009 Alexan¬ 


der et al. 20111, which select Gompton-thick AgIN at 


infrared wavelengths. Figure shows that the general 
trend is a decline of the space density of this popula¬ 
tion with decreasing redshift. Above 2 = 2, a decline 
appears towards increasing redshifts. Also, our r esults 


are in rough agreement with previous estimates (Fiore 


et al. 


et al. 


M)8l 


Alexander et al. 20111. The results of 


Fiore 


(20091 on moderate luminosity {L > 10^^-® erg/s) 


Gompton-thick AGN in the GOSMOS field is an excep¬ 
tion. The Gompton-thick AGN space density determined 
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Fig. 7. — Luminosity-dependence of the obscurer. The top row plots the obscured fraction of Compton-thin AGN (CTNOF, eq. for 
various redshift intervals. The shaded grey region are the constraints from our non-parametric method. Additionally, we compare to the 
work of|Ueda et al.|(2014|| (yellow points). F or reference, the red star symbol is placed at 70% and L = 10'^^ erg/s across the panels. In the 
top left panel, the results from local surveys ||Burlon et al.|201ll [Brightman & Nandra|2011b|| which report a similar shape. The CTNOF 
shows a distinct peak, which is placed at the red star tor focal surveys, but appears to move to sequentially higher luminosities at higher 
redshift. In the bottom row, the luminosity-dependence of the Compton-thick fraction (eq. is plotted. Our results (shaded grey) show 
that the Compton-thick fracti on is compatible with being constant at ^ 3 5% (blue star symbol for reference at L = erg/s). For 

comparison, previous surveys | |Burlon et al.|20II||Brightman Si Ueda|20I2| | are shown (see text). 




Fig. 8. — Evolution of the obscured fraction. Top panel: The 
evolution of the obscured fraction of Compton-thin AGN (CTNOF, 
eq. is plotted as a grey shaded region for L = erg/s. At 
2 ~ o, the obscured fraction was higher (75%) than t oday (50%). 
For comparison also plotted are the constraint s from |Ueda et al.| 

e (brown triangles) and |Ueda et al.|(|20I4[| (yellow points). 

n panel: The evolution ot the (Jompton-thick fraction at L = 
43.5. Our results are compatible with a cons tant Compton-thick 
fractio n of ~ 35%. We compare our results to |Brightman Sz Ueda| 
(|2012p (cyan ci rcles). Local hard X -ray Swift/BA'L observations 
(green squares, |Burlon et al.| 2011| i 
Compton-thick fraction ot ~ zU%. 


for 2 < 0.1) derived a lower 


in that study is significantly higher than our result. This 
may be attributed to contamination of obscured AGN 
samples selected in the infrared by eith er dusty star-burst 


or moderately obscured S eyferts (e.g. Georgakakis et al. 
2010 Donley et al.||201Q l. 


To study the contribution of Gompton-thick AGN to 
the total accretion luminosity output it is necessary 
to characterise the luminosity-dependence of Gompton- 
thick AGN. The bottom panel of Figure plots the 
Gompton-thick fraction, defined as 


Gompton-thick fraction := 


mH = 10 


24-26 cjn-21 


(j)[NH = 10^^“^® cm“^] 


( 3 ) 

as a function of X-ray luminosity. Within the uncer¬ 
tainties we do not find evidence for a luminosity de¬ 
pendence. If the constant-slope prior is preferred (solid 
line), a similar behaviour as in the obscured fraction is 
allowed. However, our results are consistent with a con¬ 
stant Compton-thick fraction at ~ 35%. This is the first 
time the luminosity-dependence of the Compton-thick 
fraction is constrained. 

In Figure the bottom panel shows the evolution of 
the Compton-thick fraction at L = erg/s. Our re¬ 

sults show a minor dip in the 0^1 — 2 range. This 
is due to the strong peak in Compton-thin sources (see 
above), causing a strong rise in the denominator of the 
fraction. The blue data poi nt at z < 0.1 was tak en from 
the Swift/BAT analysis of Burlon et al. (20111, where 
a considerably lower Compton-thick traction was found, 
in disagreement to our findings. This may be because 
of the small volume of our sample at those redshifts, or 
differences in the analysis and in particular the meth¬ 
ods adopted to correct for incompleteness of the AGN 
samples in the Compton-thick regime. Figure also 
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compa res our results with those of |Brightman fc Ueda| 
(20121 in the Chandra Deep Field South (cyan). Within 
the uncertainties we find broad agreement. However, 


Brightman & Ueda (20121 argued in favour of a red- 
shitt evolution (increase) of the Compton-thick fra ction 


based on their c onstraints and those of Swift/BAT (Bur- 


Ion et al. 2011). In our result we find no evidence tor 


such a trend. Our uncertainties are compatible with a 
constant Compton-thick fraction of ~ 35%. 


thin obscured, Compton-thick). One of our most strik¬ 
ing results is that Compton-thick AGN and obscured 
Compton-thin AGN contribute in equal parts to the lu¬ 
minosity density in the z = 1 — 3 range. These obscured 
AGN contribute the majority (74/^5%), while unobscured 
AGN only play a minor role. The shapes in evolution 
are fairly similar. However, the evolution of unobscured 
AGN is weaker compared to the Compton-thin obscured 
AGN when considering the z ~ 1 — 1.5 range. 



z 


Fig. 9.— The evolution of the density of Compton-thick AGN. 
Different panels correspond to X-ray luminosities brighter than 
(top), 10^®'® (middle) and 10'^'* (bottom) in units of erg/s. 
Comparing our reconstruction (black) to previous result s, we find 
overall good agreement. The result of |Fiore et al.| l [2009[ l using the 
L > 10^®'® erg/s cut is an exception. Their estimate lies higher 
compared to our results, potentially due to starburst galaxy or 
moderately obscured AGN contaminating their sample (see text). 

4.6. The contribution of Obscured AGN to the 
luminosity density 

Due to the difficulty in capturing the entire popula¬ 
tion of AGN (see previous section), the accretion his¬ 
tory of SMBHs has remained uncertain in the literature. 
Through the relation of mass accretion rate in black holes 
to luminosity, the total luminosity output of AGN can be 
used to study the total accretion history of the Universe. 
We compute the total X-ray luminosity output at any 
particular redshift, i.e. the luminosity density, by inte¬ 
grating the luminosity function as / Lx ■ <p d log Lx ■ This 
is shown in the top of Figure [To| (grey). Here, we again 
use L > erg/s, which captures essentially all of 

the luminosity density and focus on the redshift range 
z = 0.5 — 4 where our data coverage is best. The lumi¬ 
nosity density shows a broad peak in the z = 1 — 3 range, 
with a decline to both high and low redshifts. 

Additionally, in Figure we probe the contribution of 
AGN split by their obscuration (Unobscured, Gompton- 



Fig. 10.— Evolution of the X-ray luminosity density of AGN 
with Lx > erg/s, for various column densities. The lu¬ 

minosity output of AGN experiences a rise and fall in density in 
the z = 1 — 3.5 range (total as top gray shaded region). The 
strongest contribution to the luminosity density is due to ob¬ 
scured, Compton-thin (blue shaded region) and Compton-thick 
AGN (green shaded region), which contribute in equal parts to 
the luminosity. In contrast, the emission from unobscured AGN 
(red shaded region, bottom) is distinctively smaller. 

4.7. Nh distribution 

This medium luminosity interval, L « erg/s, 

where most of the evolution occurs, is also where the 
obscuration peaks according to Figure To further il¬ 
lustrate the evolution there, we look at now the column 
density distribution evolves. Figure El plots the intrin¬ 
sic Nh distribution for three redshift intervals in panels. 
Here, the boxes indicate the upper and lower 90% quan¬ 
tile on the fraction of sources in the respective Nh bin. 
The dashed line illustrates a possible Nh density dis¬ 
tribution that fits these fractions. When comparing the 
top panel to the middle panel in Figure .El. in particu¬ 
lar, the main effect appears to be that the fractions at 
Nh ~ 10^^ cm“^ increase towards higher redshifts. 

5. DISGUSSION 

We combined deep and wide-area X-ray surveys con¬ 
ducted by Ghandra and XMM-Newton to initially con¬ 
strain the space density of X-ray selected AGN as a func¬ 
tion of accretion luminosity, obscuring column density 
and redshift. Our methodology accounts for the different 
sources of errors. For example, we include in the anal¬ 
ysis uncertainties, both random and systematic, associ¬ 
ated with photometric redshift measurements. We also 
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Fig. 11. — Column density distributions at various redshifts. 
Based on our field reconstruction, we compute the fraction of 
sources in column density bins (red: “unobscured”, blue: Compton- 
thin obscured, green: Compton-thick). The boxes show the credi¬ 
ble regions containing 90% of the posterior probability. The dashed 
line illustrates a possible Njf density distribution that fits these 
fractions. 


account for the lack of redshift information for sources 
without optical or infrared counterparts. For the de¬ 
termination of the obscuration and intrinsic luminosity 
of individu al sources we use the X-ray spectral analysis 
method of Buchner et al. (20141, which takes into ac¬ 
count both the Poisson errors ot X-ray spectra and pho¬ 
tometric redshift errors. The spectral analysis also uses 
a physically meaningful, multi-component model that is 
consistent with recent ideas and observations on domi¬ 
nant AGN emission processes and the structure of matter 
in the vicinity of active supermassive black holes. An¬ 
other important feature of the analysis presented in this 
paper is the non-parametric method developed to deter¬ 
mine the X-ray luminosity function of AGN. This allows 
us to explore how the space density of AGN depends on 
luminosity, redshift, and column density without impos¬ 
ing any model. This frees us from any assumptions on 
the dependence of the AGN space density to luminosity, 
redshift, and column density. 

We find that obscured AGN {Nh > 10^^ cm“^) dom¬ 
inate the population of active super-massive black holes 


tigations (e 

g. Ueda et al. 112003 

Akylas et al.|2006 

Ueda 

et al.||20141 

Merloni &: Heinz||2013p, although our analy- 


(with uncertainties). We show that about 75% of the 
AGN space density, averaged over redshift, corresponds 
to sources with column densities Nh > 10^^ cm“^ (see 


Table . The contribution of obscured AGN to the ac¬ 
cretion density of the Universe is similarly large (75%). 
The bulk of the black hole growth across cosmic time is 
therefore taking place behind large column densities of 
gas and dust clouds. 


5.1. The role of Compton-thick AGN 

We are also able to place constraints on the number 
fraction of the most heavily obscured, Gompton-thick 
sources to the AGN population, finding it to be 38l'l7% 
of the total population. Results from previous AGN sur¬ 
veys have been divergent due to difficul ty of identifying 
Gompton-thick AG N (e.g. 15 — 20% in Akylas fc Geor- 
gantopoulos 2009 5 — 20% in Burlon et al.|20lT for the 
analyses of Swift /BAT surveys), and relatively low com¬ 
pared to the requirements from the X-ray background 
(see below). Our constraints on the Gompton-th ick frac¬ 
tion are in good ag: reement with the estimates of [Bright^ 


man & Ueda|(|2012|) of ^ 35—40%. However in their work 
they concluded an evolution of the Gompton-thick frac¬ 
tion by contrasting their high-redshift data to a local sur¬ 
vey, which reported a sig nificantly lower fraction (specif¬ 
ically Burlon et al. ]MTl| . In our analysis, we do not find 
any evidence tor a redshift evolution of the Gompton- 
thick fraction. However, our sample lacks large, shallow 
fields that can probe the local universe, and thus our re¬ 
sults may also benefit from being combined with a local 
es timate. For in s tance , the 2—10 keV X-ray based work 
of Risaliti et al. (19991 estimated that 50% of all Seyfert 
2 are observed with a Gompton-thick line of sight. If 
this estimation is taken, it is in agreement with our re¬ 
sults without the need for any evolution. However, hard 
X-ray surveys have reported signifi cantly lower v alues, 
e.g. 20j(fi% in B urlon et al. J 2 011 9 — 17% in (Bas- 


sani et al. |200D] ~ |Malizia et al.||2ul2[ [Vasudevan etalT 


2013|). This shows that estimates for a local Gompton- 
thick fraction have been diverse, making it difficult to 
make any claim of a evolutionary trend. One possible 
source of uncertainty is the sensitivity to sources with 
Nh = 10^®“^® cm“^. In this work we have explicitly as¬ 
sumed that their space density is the same as those of 
sources with Nh = 10^“^“^^ cm“^. The sparse sampling 
(only 3 secure objects in our entire sample) of these heav¬ 
ily buried population prohibits strong inferences. 

For discussing the accretion luminosity onto obscured 
AGN it is noteworthy that the fraction of Gompton-thick 
AGN does not show a luminosity-dependence (Figure]^. 
When considering currently accreting AGN, Gompton- 
thick AGN appear to play an invariant role with lumi¬ 
nosity and cosmic time. However, the large uncertainties 
in the luminosity-dependence do not allow firm conclu¬ 
sions. As Figure shows, the evolution of Gompton- 
thick AGN may follow the unobscured AGN closely in 
shape, rather than the Gompton-thin obscured AGN 
which evolve strongly. In the former case Gompton-thick 
AGN can be modelled as 35%/25% = 1.4 times more 
abundant than unobscured AGN over cosmic time and 
luminosity. In the latter case, they would be approx¬ 
imately equal in abundance and follow t he l uminosity 
dependence and its e volution (see Section 4.4 discussed 
below in Section 5.31. 


Measuring the luminosity density allows a view on the 
importance of obscuration with regard to the accretion 
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history of SMBHs. Soltan (19821 originally argued that 
a major contribution to the accretion has to be accumu¬ 
lated in obscured AGN, which were missing in surveys 
at the time. Albeit obscured AGN have been detected in 
great numbers since then, the same argument has been 
extended to Gompton-thick AGN, whose fraction has re¬ 
mained controversial. Our analysis shows the accretion 
growth of SMBHs is indeed dominated by Gompton- 
thick and Gompton-thin AGN ( 74 l 5 %), while unob¬ 
scured AGN only play a minor role. However, we also 
find that Gompton-thick and Gompton-thin AGN each 
contributing half of the luminosity output (see Figure 
10 1 , with the Gompton-thick AGN contributing 40^g%. 
This shows that Gompton-thick AGN are an important 
contributor to the accretion history. Our analysis places, 
for the first time, tight constraints on the Gompton- 
thick contribution. Results on the Gompton-thick frac¬ 
tion from studies of the cosmic X-ray background have 


varied betw een 9% (|Treister et al.|2069|) , 45%-30% (IGilli 
et al. 112007 luminosity-dependent) and 28%-60% ( (Je^ 


et al. 

2007 

et al. 

1^014 


not constrained), while |Shi et al. (2013f 
. This scatter is also due to the fact t 


re- 
that 

XRB fitting involves other parameters which are degen¬ 
erate with the Gompton-thick fra ction (e.g. additional 


reflection, see Akylas et al. 20121. It is worth empha¬ 
sizing that our estimate on the Gompton-thick fraction 
is higher than previous survey studies, but in agreement 
to those values assumed in X-ray background synthesis 
analyses. 

The overall evolution of the accretion density for the 
AGN population is presented in Figure There, the 
peak is in a broad plateau at z = 1 — 3, with a drop-off 
towards high redshift. This trend may corroborate claims 
that the accretion onto SMBHs correlates with the star - 
formation history of galaxies (see Merloni fc Heinz|2013 1. 


Unlike previous works our results do not a-priori assume 
any form of the evolution at high redshift. 

Having constrained the Gompton-thick fraction well, 
we can speculate on the origin of this obscuration. The 
fraction of Gompton-thick sources is considerably large 
(35%), necessitating that a large fraction of viewing an¬ 
gles is obscured (~ 20°, see Figure 12 1 . Gompton-thick 


obscuration is associated with the torus, as such column 
densities are not typically reached by galactic gas. As 
a G ompton-thick, smooth obs curation would be unsta¬ 
ble (Krolik & Begelman 19881, the current working hy- 
pothiesis is that the obscuration comes in discrete clouds. 
Under this clumpy torus model obscured views are pro¬ 
duced when clouds are encountered in the line of sight. 
If the Gompton-thin obscuration is also associated with 
the torus, Gompton-thick AGN are a simple extension of 
Gompton-thin obscured AGN. In principle it would also 
be possible then that Gompton-thick AGN have just a 
larger number of clouds in a clumpy torus. In this case, 
their cumulative line-of-sight obscuration would provide 
a Gompton-thick view. Alternatively, these objects may 
have denser clouds overall. Another possibility is that the 
cloud density is distributed such that both regimes are 
covered. For further research in this regard, the unbiased 
column density distributions shown in Figure [^can pro¬ 
vide observational constraints for clumpy torus models. 
The alternative scenario is that Gompton-thin obscura¬ 
tion is associated with a medium of different extent than 


Gompton-th ick obscura tion (e.g. through galactic and 
nuclear gas, Matt| 20001. We discuss both kind of mod¬ 
els in Section~ f5.3[ '\ 
dependence. 


when also considering the luminosity 


5.2. Obscuration-dependent evolution 

In Figure we show that the space density of ob¬ 
scured Gompton-thin AGN experience a much stronger 
rise in the redshift interval z = 0.5 — 1.2 compared 
to unobscured AGN. This demonstrates an observed 
obscuration-dependent evolution. Our method explicitly 
contains a smoothness prior, preferring that the space 
density of AGN as a function of column densities is the 
same. In contrast, the results show strong differences in 
the evolution of unobscured and obscured Gompton-thin 
AGN. We can thus conclude that this result was driven 
by the data. 

In Figure we find that the fastest evolving popu¬ 
lation is the one with column densities Nh = 10 ^^ — 
10^^'^ cm“^. A possible interpretation of these trends is 
that different levels of obscuration correspond to media 
with different spatial extents. We investigated this fur¬ 
ther by focusing on the behaviour of moderately obscured 
AGN {Nh = 10^2 - lO^^ cm-^). 

5.3. Luminosity-dependence and evolution of the 
obscured fraetion 

An important result from our analysis is that the ob¬ 
scured Gompton-thin AGN fraction (GTNOF, eq. de¬ 
pends on both luminosity and redshift (see Figure [^. 
The luminosity dependence can be described by a peak of 
the GTNOF at a given luminosity followed by a decline at 
both brighter and fainter luminosities. The redshift de¬ 
pendence is manifested by a shift to brighter luminosities 
of the peak of the obscured AGN fraction with increasing 
redshift. The shift of this peak causes observers consid¬ 
ering mostly AGN with luminosity L > 10“^^ erg/s to see 
an increase in the fraction of obscured AGN with redshift 
(as shown in Figure]^. Previous studies also find that 
the fraction of obscured AGN decreases with increasing 
luminosity (e.g. Ueda et al.|2003 Akylas et al.|2006 Ueda 
et al.|20f4 |. When we consider a luminosity range where 
we have good constraints {L = 1043 - 2-46 0 j.g/gj^ ij^e ob¬ 
scured fraction rises with redshift so that at z > 2.25, 
83^1% of AGN are obscured {Nh > lO^^ cm”^)^ as com¬ 
pared to the local z = 0 value z = 75 ^ 4 %. 

Our analysis also establishes a decline of the obscured 
AGN fraction at low luminosities. This trend has only 


Bur Ion et al. 2011 
also leads to the 
consider only in- 


been found in the local Un iverse dBm 
Brightman fc Nandra|| 20 lTb l. Figured 
conclusion that observers who wouldc 
trinsically faint AGN {L < 10'*^ erg/s) would observe a 
decrease in the fraction of obscured AGN with redshift. 
In other words, the magnitude of the trend is determined 
by the sample selection. 

Obscurer models that do not invoke any obscuration- 
dependence are ruled out by our findings. These include 
the simple torus where the obscuration distribution is 
produced purely geometrically. Despite these shortcom¬ 
ings, this picture may serve as a visualisation of the ob¬ 
served trends. The luminosity-dependence and evolution 
of obscured/Gompton-thick AGN is illustrated in Figure 
12 If we assume a toroidal geometry for the obscurer 
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low-L mid-L high-L 


f low-z high-z 

U 

Fig. 12. — Illustration of the population-averaged obscuring 
toroid. The obscured and Compton-thick fraction and their lu¬ 
minosity dependence is shown using a corresponding opening an¬ 
gle above which the sky is obscured to the X-ray emitter (center). 
The luminosity dependence is depicted for low-redshift (left) and 
high-redshift (right) AGN using radial shells: as the luminosity in¬ 
creases, more sky is visible. Notice that the difference between a 
low-redshift and high-redshift source does not affect the lowest and 
highest luminosities, but only where the transition occurs (blue 
arrow). The opening angle corresponding to the Compton-thick 
transition remains roughly constant. For accurate numbers, refer 
to Figures [^and|^ 


shared by the whole population, we can define an open¬ 
ing angle that reproduces the obscured/Compton-thick 
fraction. The luminosity dependent increase of this open¬ 
ing angle from 30° to 40° is depicted using radial shells 
for low-redshift (left) and high-redshift (right) AGN. No¬ 
tice that the difference between a low-redshift and high- 
redshift source does not affect the lowest and highest 
luminosities, but only where the transition occurs (blue 
arrow). Thus, it is not the case that high-redshift sources 
merely have more obscuration overall. If this were the 
case, a boost in the obscured fraction would be visible 
over all luminosities. Instead the peak of the obscured 
fraction has moved to lower luminosities, thus reducing 
the obscurer in moderately bright AGN (see Figure [^for 
accurate numbers). 

We now consider a few exemplary models that can 
qualitatively account for the decrease of the obscured 
fraction towards low luminosity or the decrease towards 
high luminosities. An additional uncertainty is whether 
the relevant obscuration occurs on small or large scales, 
i.e. in the torus or on galactic scales. In the follow¬ 
ing, we present for both scenarios a model each for the 
high-luminosityand the low-luminosity effect. These are 
illustrated in Figure [T^ 

The receding torus (RT) model is often invoked to ex¬ 
plain the observed decrease of the obscured AGN fraction 
at bright luminosities. This scenario postulates a toroidal 
geometry for the obscurer with opening angle that in¬ 
creases with increasing luminosity. This is attributed 
to photon pressure pushing away the obscuring mate¬ 
rial, photo-ionisation of the gas clouds or sublimation 
of the dust by the photon field of the active black hole 


dLawren^ 1991 


, Simpson|2005| 

los|2008 ). such a variation ot the scale height ot the torus 

has been claimed via infrared observations dMaioli no fc 
Rieke 1995, Lusso et al. 2013, Toba et al.J 20l4] r In 


Akylas k. Georgantopou- 


the c ontext of clumpy to rus models (H5nig fc BeckerT] 


2007| |Nenkova et alT2008|), a more intense radiation field 
al 


leads to the (partial) ionisation of individual dust and 
gas clouds thereby increasing the effective solid angle of 
unobscured sight-lines. This model (RT) thus represents 
a direct, causal connection between the X-ray luminosity 
and the unobscured line of sight. This scenario is illus¬ 
trated in the top right of Figure In the RT model it 
is, however, difficult to explain the observed evolution of 
the luminosity dependence. The main evolutionary effect 
observed is an increase of the turnover luminosity with 
redshift, causing the onset of the high-luminosity effect 
at subsequently higher luminosities (see Figure Q. For 
a purely nuclear obscurer scenario such as RT, no such 
evolutionary effect is predicted. 

Alternatively, it is also possible to interpret the ob¬ 
served luminosity dependence of the obscured AGN frac¬ 
tion in the context of evolutionary models, in which dif¬ 
ferent levels of obscuration roughly correspond to dif¬ 
ferent stages of the growth of supermassive black holes. 


of gas inflows caused by mergers (e.g. 

Hopkins et al. 

2006a 

Somerville et al. 

[20081 or secular processes (e.g. 

h'anida 

Ids et al. 1120121 |J 

Sournaud et al.j 

20071 [Giotti &[ 

Ustriker 

19971 L 

during this period super-massive 

black ho 

es grow tast by accreting material close to the 


Eddington limit and are also typically obscured by the 
inflowing material. Eventually however, the energy out¬ 
put of the central engine becomes sufficiently powerful 
to drive outflows, which can blow away the obscuring 
clouds of dust and gas. The central engine then shines 
unobscured for a brief period before its luminosity out¬ 
put declines as a result of the depletion of the available 
gas re servoirs. It is further proposed (e.g. [Hopkins et ^ 


2005 2006b) that luminous AGN brighter than the break 


of the X-ray luminosity function are dominated by fast 
accreting black holes close to the peak of their growth 
phase. Lower luminosity sources include a large fraction 
of AGN during the decline stage of their activity. 

In the above scenario the decrease of the obscured AGN 
fraction with increasing luminosity may be linked to the 
blow-out stage of AGN. This scenario (EB) illustrated in 
the bottom right panel of Figure [^ At brighter accre¬ 
tion luminosities it is more likely thaf AGN-related feed¬ 
back mechanisms become more efficient, thereby sheding 
the gas and dust cocoons in the vicinity of supermas¬ 
sive black holes. Lower accretion luminosities, below the 
break of the X-ray luminosity function, include a large 
fraction of AGN past the peak of their activity, when 
the gas reservoirs that feed and obscure the black hole 
have already been depleted. Therefore, one might expect 
a larger fraction of unobscured AGN at low accretion 
luminosities. The evolutionary picture outlined above 
is therefore consistent with the observed luminosity de¬ 
pendence of the obscured AGN fraction: a maximum at 
a given luminosity and a decline at both brighter and 
fainter luminosities. 

In the EB model, the observed redshift dependence of 
the turnover luminosity is also difficult to explain. Un¬ 
less the physical conditions of the obscuring material are 
a function of cosmic time, it is not obvious how the phys¬ 
ical mechanisms that affect the geometry of the obscurer 
on small-scales can produce a luminosity dependent ob¬ 
scured AGN fraction that changes with redshift. 
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Fig. 13.— Four models that address the luminosity-dependence of the obscured fraction. This categorisation distinguishes by the 
luminosity-dependent effect, the increase of the obscured fraction (left column) and the decrease of the obscured fraction (right column). 
The top of Figurej^is repeated here and simplified. In the four models considered (see text for details), the relevant obscurer is either the 
torus (top row) or galactic obscuration (bottom row). In the illustration of the models, the red star represents the X-ray source, while grey 
clumps illustrate the cold obscurer. 


A rescuing argument for the EB scenario is possible. 
It is now fairly established that the black hole mass of 
the a verage, currently accreting AGN incre ases with red- 
shift ( MerlonipOOd Merloni fc Heinz|2013 1. This means 


that at high redshift, the typical black hole undergoing 
accretion has already accumulated a large mass, while 
at low redshift, small black holes undergo accretion. Re¬ 
cent simulations have shown that such a anti-hierarchical 


mation fe.g. IFanidakis et al. 

2012 

Hirschmann et al. 

2012 

Enoki et al.|2U14p, and can exp 

lain the decrease of 


nario, the peak of the luminosity function is linked to the 
Eddington luminosity, as brighter AGN remove their ob¬ 
scuration. If the black hole mass increases with redshift, 
so does the Eddington luminosity, causing the on-set of 
the high-luminosity decrease of the obscured fraction to 
occur at higher luminosities. The EB scenario is thus 
compatible with our observations, if the average black 
hole mass increases with redshift. 

This model can affect both the small-scale, nuclear ob¬ 
scurer as well as the gas in the galactic vicinity of the 
black hole, and our observations may be incapable of 
distinguishing the effects. For this reason, the illustra¬ 
tion of this effect (EB) in Figure 13 might be placed in 
both categories. 

Lets now consider the low-luminosity effect, which is a 
decrease of the obscured fraction towards low luminosi¬ 
ties. The evolutionary scenario might also be capable of 
explaining the low lum inosity increase of the obscuration 
(Hopkins et al. 2006aI. If the observed column density 
is due to galactic streams, feeding the black hole simul¬ 
taneously obscures it. One would then expect a simple 
correlation of the accretion luminosity and the obscura¬ 


tion, at least at luminosities below the Eddington limit. 
This effect (Gas Available, GA) is illustrated in the bot¬ 
tom left of Figure [T^ 

Under the GA scenario (bottom left of Figure [I^, the 
accretion luminosity is simply correlated to the available 
gas. While observations have shown an overall reduction 
of the gas and dust c ontent of galaxies over cosmic time 
(Tacconi et al. 20131, under the GA model this would 
increase the total number of AGN at high luminosities, 
but the luminosity-dependence of the obscured fraction 
would remain unaffected. In particular, the GA model 
does not predict the observed increase of the turn-over 
luminosity as suggested by Figure 0 

The low-luminosity increase of the obsc ured fraction 
has been observed before in local samples dBurlon et al. 


2011 Brightman & Nandra 2011b Elitzur &: Ho 2009|. 
hbr these sources, it is more probable that the obscurer 
is nuclear, as many of these sources provide the opportu¬ 
nity to study both the galaxy and the infrared emission 
of the torus. One model is that the obscurer is caused 


by clumps in the disk wind (Elitzur & Shlosman 20061. 
This scenario was originally created to explain the verti¬ 
cal supp ort for a cold dusty torus t hat would otherwise 
collapse (Krolik & Begelman 19881. In the disk wind 
model, the accretion disk is only capable of projecting 
clouds beyond a certain luminosity, predicting the ab¬ 
sence of a clumpy torus in low-luminosity AGN. This 
scenario of cloud production (GP) is illustrated in the 
top left of Figure [l^ 

While the cloud production model reproduces the nec¬ 
essary increase of the obscured fraction, it does not pre¬ 
dict any evolution over cosmic time when taken at face 
value. However, under this model, the critical luminosity 
below which the obscurer can not be sustained is strictly 
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determined by the black hole mass (Elitzur & Shlosman 
200^|Elitzur|®l. This model has received observa- 


h- 

Eli 


tional support in |Elitzur & Ho (20091 who observe a 
black hole mass dependence ot the luminosity depen¬ 
dence from a local survey. This model, taken together 
with the black hole mass evolution explained above, can 
thus explain the observed evolution. 


TABLE 4 

Predictions of the discussed models. 


ModeF 

Obscuration*^ 

Low- 

L 

effect^ 

High- 

L 

effect‘d 

Evolution® 

RT 

nuclear 

/ 


BHMf 

CP 

nuclear 


/ 

no 

EB 

nuclear or galactic 


/ 

BHMf 

GA 

nuclear 

/ 


no 

RT+CP 

nuclear 

/ 

/ 

not for high-L 

GA+EB 

galactic 

/ 

/ 

not for low-L 

CP+EB 

nuclear or both 

/ 

/ 

BHMf 


^ Model name as presented in Figure [Tsl 
Whether the luminosity-dependent layer of Obscuration is 
associated with galaxy-scale obscuration (“galactic”) or the 
“torus” (“nuclear”). 

“ Whether the model predicts a decrease of the fraction of ob¬ 
scured AGN towards low luminosities. 

Whether the model predicts a decrease of the fraction of ob¬ 
scured AGN towards high luminosities. 

® Whether the model predicts any evolution of the effect with 
cosmic time. 

f The model scales with black hole mass. Evolution is predicted 
if the average accreting black hole mass changes over cosmic 
time. 


To conclude, we have considered four effects (simple 
models) which are summarized in Table Individually, 
they only partially explain the observations. However, 
the CP model, taken together with the EB scenario can 
explain the observations if the black hole mass function 
evolves with redshift. Under this view, the obscurer is 
vertically extended above a certain luminosity, and the 
torus disappears when the Eddington luminosity drives 
away the obscurer, unveiling the bright X-ray source 
for a (cosmically) brief time. Here it remains uncertain 
whether the obscuration removed is galactic, nuclear, or 
both. Both effects engage at a luminosity that is depen¬ 
dent on the black hole mass, which provides the cosmic 
evolution effect. 

Future observations may shed light on the dependence 
of the galactic gas with the AGN X-ray luminosity. 
For instance, AL MA (Atacama Larg e Millimeter/Sub- 
millimeter Array, Beasley et al. 20061 measurements of 
gas fractions in host galaxies ot unobscured and obscured 
AGN can be compared. If they do not differ, the GA 
model is ruled out. A similar approach by comparing lu¬ 
minous and faint AGN can also probe the predictions of 
the EB model, and establish whether the obscuration is 
nuclear. Under the EB model, gas motions are expected 
in bright AGN, which may be detectable wit h future 


X-ray spec troscopes on the ATHENA mission (Nandra 


et al. 20131. Furthermore, to establish that the evolution¬ 


ary trend is caused by black hole mass evolution, black 
hole mass-matched samples of luminous AGN may be 


considered. If such a sample shows no redshift evolution 
in the luminosity-dependence of the obscured fraction, 
then the trend is due to black hole mass differences. The 
GANDELS fields with their high-resolution near-infrared 
imaging is the best candidate for such a study. 

The models as presented here are a rudimentary de¬ 
scription of the involved physical processes and re¬ 
quire further refinement in their predictions and self- 
consistency, e.g. via numerical simulations. Among the 
challenges in this regard is that the nuclear ob scuration 


can not be resolved in evolutionary models (e.g. Hopkins 
et al.|2006a, Somerville et al.|2008 1 and thus the torus is 


not yet treated selt-consistently. Recently, it has become 
possible to self-consistently model the radiative a nd hy¬ 
drody namic processes that maintain the torus (e.g.|Wada 


20121. Further research is needed in terms of whether the 


torus can, by itself, reproduce the luminosity dependence 
of the obscured fraction. Also, it has to be clarified which 
processes maintain a cold, Gompton-thick torus, with the 
geometric extent implied by a 35% Gompton-thick frac¬ 
tion. 


6. GONGLUSIONS 

We combined deep and shallow, wide-area X-ray sur¬ 
veys conducted by Ghandra and XMM-Newton to con¬ 
strain the space density of X-ray selected AGN as a func¬ 
tion of accretion luminosity, obscuring column density 
and redshift. An important feature of our analysis is 
the non-parametric method developed, which does not 
require any assumptions on the shape of the luminosity 
function or its evolution, allowing the data to drive our 
results. Furthermore, we take into account all sources of 
uncertainties, allowing us to robustly constrain the evo¬ 
lution of unobscured and obscured AGN, including the 
most heavily obscured Gompton-thick. 

We find that obscured AGN with Nh > 10^^ cm“^ 
account for 77^g% of the number density and 74)'(g% of 
the luminosity density of the accretion SMBH popula¬ 
tion averaged over cosmic time. Gompton-thick objects, 
with Nh > 10^"^ cm“^ account for approximately half 
the number and luminosity density of the obscured pop¬ 
ulation, and 38)'l7% of the total. 

There is evidence that the space density of obscured, 
Gompton-thin AGN evolves stronger than the unob¬ 
scured or Gompton-thick AGN. This is connected to 
the luminosity-dependent fraction of obscured AGN. At 
higher luminosities, fewer AGN are obscured. However, 
at higher redshift, this effect sets on at significantly 
higher luminosities. In the luminosity range used in 
this study, the fraction of obscured AGN increases from 
75^4% to 83)'l3% at z > 2.25. This is due to only a small 
luminosity range around L* w 10^"^ erg/s which changes 
its obscured fraction. In particular, the space density 
evolves fastest around Nh ~ 10^^ cm“^. In contrast the 
fraction of Gompton-thick AGN relative to the total pop¬ 
ulation is consistent with being constant at ~ 35% inde¬ 
pendent of redshift and accretion luminosity. The contri¬ 
bution to the luminosity density by Gompton-thick AGN 
is 40)'lg%. The robust determination of a large fraction 
of Gompton-thick AGN consolidates AGN X-ray surveys 
and studies of the cosmic X-ray background. It also im¬ 
plies a physically extended nuclear obscurer. 

We also find that the fraction of moderately obscured 
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AGN decreases not only to high but also to low lumi¬ 
nosities, in qualitative agreement with findings of lo¬ 
cal surveys. We have discussed the observed trends 
of the obscured fraction with Lx and z in the con¬ 
text of models that either assign obscuration to the 
torus or galaxy-AGN co-evolution effects, with varying 
luminosity-dependent effects. Both classes of models can 
qualitatively explain our results but require that SMBH 
evolve in a downsizing manner, i.e. larger black holes 
form earlier in the Universe. 
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APPENDIX 

A. LUMINOSITY FUNCTION ANALYSIS 

We here go through a derivation of the statistical foot¬ 
ing of analysing population d emographics b y revi ewing 


and combin ing the works of Loredo 


(20041 and Kelly 


et al. (2008). Their main difference is whether the Bino¬ 


mial distribution or its approximation, the Poisson dis¬ 
tribution, is used. 

We intend to estimate the number density per comov¬ 
ing volume (Mpc^) of AGN, as a function of various prop¬ 
erties, specifically X-ray luminosity, redshift and obscur¬ 
ing column density. This extended luminosity function 
will describe the evolution of the X-ray population and 
all sub-populations (e.g. Compton-thick AGN). The dif¬ 
ficulty is that each of these properties influences our abil¬ 
ity to detect objects with such properties. Lets assume 
we can describe the probability to detect a source with 
properties C as p{'D,C). 

After we analysed each object in detail, we can bin our 
sample in so small bins that only 1 item at most can be 
in each bin Ci. Then just applying the Binomial/Poisson 
distribution, assuming that in k bins a source is detected 
(and no detections in the other n — k bins), yields: 


B{k]n,p) = xY[p{V,d,,C,)x (Al) 


]Jp(X>,fii,C*) 


i—k 


P(fc;n,p) = — X JJn X p(T>,dj,Cj) X (A2) 


exp < -y^p{'D,di,Ci) 


i=l 


The first equation gives the likelihood based on the 
Binomial distributions, and here p{'D,di,Ci) denotes the 
probability of not detecting a source with properties Ci 
and having obtained the observed data di. In the Poisson 
formula, only p{'D,di,Ci) occurs, which instead denotes 
the probability of detection. 
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W e no tice that the term marked with a star in Equa¬ 
tion |A2| remains the same, regardless of the number of 
detections k in the sample (only dependent on n). The 
sum in this term therefore has to be always the same, 
independent of the specific sample di. In the Poisson 
formalism, the exponent of this term has a specific mean¬ 
ing - it refers to the expectation value of the sample (a 
a-priori assumption). We can thus remove the depen¬ 
dency of the data and treat it as an a-priori probability. 
Mathematically, we integrate over all possible data di. 


P(fc;n,p)=exp X (A3) 

k 

Y\{nxp{Vi,di,Ci)) (A4) 

We used YJi=iP{'^i^di,Ci) = J2'^^iP{T>i,Ci). By 
analogy, for the binomial distribution, we use 
Jp{'Di,di,Ci) ddi = p(Vi,Ci). Mathematically, we again 
integrate over all p ossible data dj , b ut for all non- 
detected sources (see Kelly et al. (20081). 


B{k;n,p) = x ]Jp(X>j,C*) x (A5) 

i=k 
k 

l[p{V„d„C,) (A6) 

= (”) X exp |^lnp(I?,,Ci)+ (A7) 

InpiVi.di.Ci) 

i=l 

Now we replace our n discrete bins by a continuum. 
The probability p{'Di,di,Ci) can then be non-zero over a 
range of the parameter space. After all, we do not know 
the true parameter C. Thus, we replace '^^p{Vi,di,Ci) 
by jgp{Vi,di\C)dC. 




P{k-, n,p) = — X exp j P(^iC) dC+ 

(A9) 

k ^ 

y^ln [p{'Di,di,C)dc\ 

) 

(AlO) 

B{k;n,p)= x exp{ 

(All) 

(n — fc) • In J p{'D, C) dC + 

(A12) 

y^ln J piV^,di,C)dC 1 

(A13) 


We will now expand p{'D, C) and p{T>i, di,C) using con¬ 
ditional probabilities and discuss the meaning of the oc- 
curing terms. 


p{V,C)=p{C)-p{V\C) (A14) 

p(T>,di,C)=p{C) -pidilC) ■p{'D\di,C) (A15) 


Here, p{C) is the probability of finding an object with 
characteristics C {{log L, z, log Nh} here). This is essen¬ 
tially the luminosity func tion we are interested in. The 
term p{'D\C) in Equation A14 denotes the sensitivity to 
such objects. 

In Equation |A15 the second term, p{di\C), is related 
to the spectral analysis. It denotes the likelihood that 
this data was generated from a source with e.g. lumi¬ 
nosity L and column density N^. Finally, we have to 
consider the probability of detecting this object, given 
the data and characteristics C. This is almost surely 1, 
as having data associated with an object implies having 
detected it. However, there is a subtlety here that has 
been overlooked so far. Some astronomers use the sen¬ 
sitivity p{T>\C) ~ specifically the area curve of where the 
obje ct was detecte d - here, in place of p{'D\di,C) = 1. 
But Loredo| (20041 makes a strong point arguing that 
the sensitivity function must not be used here. Why is 
it done in practise anyways? 

The X-ray source detection algorithms typically use a 
area of a certain, relatively small size, which contains 
70% of the energy contributed by the source. The back¬ 
ground for this region is estimated, and the probability 
for the background to produce the observed number of 
counts computed (no-source probability). If a threshold 
is exceeded, the spectrum is extracted. In this step a 
much larger region is used, containing more background 
counts. While in the small detection area it was not pos¬ 
sible for the background to produce the observed counts, 
for faint sources it may be possible for the background 
to produce the observed counts in the spectrum extrac¬ 
tion region. This is because the background contribution 
grows linearly with area, but the source contribution al¬ 
most stays the same. The likelihood of our analysis thus 
does allow having no flux from a faint source, i.e. lumi¬ 
nosity zero. This contradiction to the detection proba¬ 
bility stems from the fact that we have not used the in¬ 
formation that the counts are concentrated in detection 
region. In this case, p{V\di,C) ^ 1. 

Specifically, we could write it as p{'D\di,C) = p{> 
k\di,C), the probability of having more than k counts in 
the detection region out of the extracted spectrum counts 
di, where k denotes the number necessary for detection 
at the current position. An approximation to this num¬ 
ber is p{> k\C), the probability to produce the number 
of counts required for a detection at the current position. 
A further approximation is f p{> k\C) dA/A = p{'D\C), 
which is the area-average sensitivity curve. Given this, 
we understand now why some works use the area curve 
in the data term, even though it should not be done: it 
is an attempt in fixing a loss of information introduced 
when the detection and data analysis processes differ. 

We now put all the information together and arrive at 
the relevant likelihoods: 


B{k; n,p) 



p{V\C)-p{C)dC] >(A16) 


%—k 
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k p 

n / p{C)-p{d,\C)-p{V\d,,C))dC 


P{k;n,p) = — X exp 


piC) ■ p{V\C) dC } >(A17) 


k p 

n / p{C)-p{d,\C)-p[V\d,,C))dC 


Comparing Equation 9 in Loredo 


^ ^ (|2004| with Equa¬ 

tion mi above, we have obtained tire same result with 
the notation key p{di\0) = li, p{'D\di,C) = 1, C = m, 
n = LoSm and p(I?|C) = ri(m). The derivation is similar, 
however we started from first principles and continued in 
a sound and complete way for both the Poisson and the 
Binomial distribution. 

Comparing Equation 5 in Kelly et al.| (20081 with 
Equation |A16| above, we have obtained the same result 
with the notation key C = {Lj,Zj), p{'D\C) = p{Ij = 
0\Lj,Zj), p{C) = p(£,-,2:d^)i p{T>\di ,C) = 1 as well as 
p{di\C) = 1, since Kelly et al.f (|2008|) neglects any mea¬ 
surements uncertainties in the derivation, assuming that 
luminosity and redshift can be determined perfectly. 

Finally, we summarise the logarithms of the likeli¬ 
hoods, neglecting constants. 


In/is = (n — fc) X In Jp(T>\C) ■ p(C) dC + (A18) 

k p 

^In / p[C)-p{d,\C)-p{V\d,,C)dC 


In /ip = - / p{C) ■ p{V\C) dC + (A19) 


^In J piC)-pid,\C)-p(V\d,,C)dC 

i—1 


The gene ral-purpose f ramew ork presented here is ad¬ 
vocated by Kelly et al. (2008). It is a general approach 
for the inference ot population demographics based on 
specific samples. The importance of this step - rather 
than analysing trends based on sample statistics - and 
incorporating selection biases, can not be over-stated. 

In this work, we restrict ourselves to the Poisson like¬ 
lihood. This allows us to use the space density directly 
(rather than separating sampled volume and probabil¬ 
ity). For p(C) we insert the luminosity function model 
(j){logL, z, logNn) X For p(I?|C) we insert the area 
curve. For p{'D\di,C) we use the area curve of the spe¬ 
cific field as an approximation (see above). For p(di\C) we 
should use the likelihood of the spectral analysis. How¬ 
ever, we did use intermediate priors in the data analy¬ 
sis. We thus extend p{di\C) = 


term is the posterior distribution computed in the spec¬ 
tral analysis, and the prior has to be divided away again. 
As we used flat priors in logL, log Nh, and z (subject to 
additional information from other wavelengths), which 
are the units of the integral over C = {logT, z, log A'//}, 
this division only contributes to the likelihood as a fixed 


offset, and is not relevant for further analysis. If different 
intermediate priors had been used, here they would be 
have to be divided away. In fact, we employed a Gaussian 
prior on the photon index T. The luminosity functions (j) 
we consider however are in fact (/<(log L, z, log Nh) x /(T) 
where / is the Gaussian prior used, thus we also do not 
require a correction here. 

Regarding the practical evaluation of the likelihood 
function, we use posterior samples from our spectral 
analysis as just described. For the first integral, which 
describes the expected number of detected sources, we 
employ a fixed grid, as these are fast to compute but 
also ensure a consistent estimate between likelihood eval¬ 
uations. We use 40 grid points in each dimension. By 
pre-computation of all weights except for the luminos¬ 
ity function (j), the likelihood function reduces to (fast) 
addition and multiplication operations. 


B. ROBUST PHOTOMETRIC REDSHIFT 
PROBABILITY DISTRIBUTIONS 

In this section we elaborate on incorporating system¬ 
atic uncertainties into photometric probability distribu¬ 
tions from SED fitting. 


The photometric fitting me thod of Salvat o et al. ( 2011 


2009), using software from Ilbert et al. 




IS capable of not just computing a most likely redshitf 
point estimate, but producing red shift probability distri - 
butions. Previous studies, such as Buchner et al. (2014), 
used the photo-z probability distributions directly to in¬ 
corporate the uncertainty of the redshift estimate. How¬ 
ever, additional to the uncertainty due to measurement 
errors, the method of photometric redshift has system¬ 
atic errors due to incomplete template libraries, failures 
in automatically extracting correct fluxes at the edge of 
images, blending of sources causing the flux to be convo¬ 
luted, etc. The most important systematic contribution 
is probably incorrect associations. 

On a high-level, the quality of redshift point estimates 
for a specific sample can be described by two quantities, 
which are estimated from a sub-sample where spectro¬ 
scopic redshift s are available: the outlier fraction ry and 
the scatter cr (Salvato et al. 2011). The outlier fraction 
describes the failure rate or accuracy of the redshifts, i.e. 
how frequent the redshift estimate is off by far: 


r] : fraction where 


•^phot 


1 


> 0.15 


The scatter width cr describes how precise the redshifts 
are scattered around the true value. The sample can be 
characterised using the normalised median absolute devi¬ 
ation (NMAD, Hoaglin et al. 19831 as a robust estimator: 


o’NMAD = 1-48 X median 


■^■phot 


1 -|- Z, 


spec 


These two quantities are interpreted that the redshift 
estimates Zphot are distributed along a normal distribu¬ 
tion centred at the true value Zgpec with standard devi¬ 
ation (Tnmad ■ (1 + ^spec)- Additionally, a fraction ry of 
the redshift estimates are drawn from a different distri¬ 
bution, which can be described as flat over a redshift 
range (0, Zmax)- Here we considered z^ax = 5 for out¬ 
liers. The probability distribution of the point estimator 




















































23 


1 data ■( 

2 int indices [27344,3] ; 

3 real weights [27344] ; 

4 vector[ll-l] widthsO; 

5 vector[14-l] widthsl; 

6 vector[7-l] widths2; 

7 int lengths [2046] ; 

8 int chain_indices [780,3] ; 

9 real chain_weights [780] ; 

10 } 

11 parameters { 

12 real<lower=-40,upper=0> y[ll,14,7]; 

13 > 

14 model { 

15 real sigmaL; 

16 real sigmaz; 

17 real sigmaNH; 

18 real dataterms[2046]; 

19 real detectionterms [780] ; 

20 real loglike; 

21 

22 sigmaL <- 0.5; 

23 sigmaz <- 0.5; 

24 SigmaNH <- 0.75; 

25 

26 /* L smoothness prior: 2nd derivative is small */ 

27 f or (i in 3 : 11) { 

28 for (j in 1:14) { 

29 f or (k in 1 : 7) { 

30 y[i,j,k] ~ normal((y[i-1,j,k] - y[i-2,j,k]) * widths0[i-l] / widths0[i-2] + y[i-l,j,k], sigmaL); 

31 } } } 

32 

33 /* z smoothness prior: 2nd derivative is small */ 

34 f or (i in 1:11) { 

35 for (j in 3:14) { 

36 f or (k in 1 : 7) { 

37 y[i,j,k] ~ normal ((y [i , j-1 , k] - y[i,j-2,k]) * widthsl[j-l] / widthsl[j-2] + y[i,j-l,k], sigmaz); 

38 } } } 

39 

40 /* NH smoothness prior: 1st derivative is small ♦/ 

41 f or (i in 1 : 11) { 

42 for (j in 1:14) { 

43 f or (k in 2 : 7) { 

44 y[i,j,k]'' normal (y[i ,j ,k-l] , sigmaNH) ; 

45 } } } 

46 

47 { /* individual objects */ 

48 int m; 

49 m < - 0; 

50 for (i in 1:2046) { 

51 int a; int b; int c; 

52 real v [lengths [i]] ; 

53 

54 for (j in 1:lengths [i]) { 

55 int 1 ; 

56 1 < - m + j ; 

57 /* look up interpolation point and compute weights */ 

58 a <- indices[l,l] + 1; 

59 b <- indices[l,2] + 1; 

60 c <- indices[l,3] + 1; 

61 v[j] <- exp(y [a,b,c]) * weights[l]; 

62 } 

63 m <- m + lengths [i] ; 

64 /* add value to likelihood */ 

65 dataterms[i] <- log(sum(v) / lengths[i]); 

66 } 

67 } 

68 

69 /* detection integral */ 

70 for (k in 1:780) { 

71 int a; int b; int c; 

72 real v; 

73 /* look up interpolation point and compute weights */ 

74 a <- chain.indices[k,1] + 1; 

75 b <- chain.indices [k,2] + 1; 

76 c <- chain.indices[k,3] + 1; 

77 V <- chain.weights[k] * exp(y [a,b,c]) ; 

78 detectionterms [k] <- v; 

79 } 

80 loglike <- sum(dataterms) + sum(detectionterms) ; 

81 increment.log.prob(loglike); 

82 } 

Algorithm 1: Stan code for estimating the field. The contant-slope prior implementation is shown. We prepared the 
integrals to be computed using importance sampling chain points evaluated the respective field bin value (indices) 
and multiplied by the weights. 
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incorporating systematic errors (SYSPE) can be written 
as 


SYSPE(z) := ry X U(0, + (1 - 7 ?) x N(zphot, a) 

This modelling in fact matches the empirical distribution 
very well when considering the cumulative distributions 

of l-Sphot •^spec 1/(1 + '^(spec) ■ 

In this work, we would like to similarly incorporate 
catastrophic outliers. These are not contained in the 
photometric redshift distribution (PDZ) from SED fit¬ 
ting. Thus, we follow an analoguous modeling of the 
uncertainty contribution: We allow a broadening of the 
PDZ by convolution with a Gaussian kernel of width a 
and add a flat probability plateau with weight fj: 

SYSPDZ(z) :=^xU(0, z„,ax) + (l-f 7 ) x (PDZ * N(0, a)) 

Note that for the special case a — 0 and fj — 0, SYSPDZ 
becomes PDZ. The two parameters a and fj have a 
slightly different definition than a and 77 above, but take 
the same roles for characterise systematic uncertainty. 
We fit for the parameters using the spectroscopic sub¬ 
sample by maximising the likelihood 

£((7,77) =[] SYSPDZ,(z,pec..) 


defined as the product of the value of the modified PDZ 
at the true redshifts ^spec- PDZs that are far off the true 
value will demand an increase in a and fj, however raising 
both values diminishes the SYSPDZ value overall due to 
normalisation of probability distributions. 

We And the best fit parameters a = 0.024, 0.048, 0.029 
and fj = 1.8%, 0.6%, 2.3% for the CDFS, AEGIS-XD 
and G-GOSMOS samples respectively. These numbers 
are not directly equivalent to t he point estimate b ased 
definition of gN MAn and 77 (see Salvato et al. 2011 and 


Hsu et al. 2014 1 , but the numbers are comparable. We 
thus use tne smoothed SYSPDZ distributions instead of 
the PDZs. 

A simpler method for incorporating the systematic un¬ 
certainty would be to use the modelling based on the 
point estimator (SYSPE). But when comparing the like¬ 
lihood value, we And that SYSPE always has lower like¬ 
lihood values than the kernel smoothing SYSPDZ. This 
shows that the PDZ contains valuable information due 
to e.g. secondary solutions, which is missing in the point 
estimator Zphotj and that SYSPDZ is the most faithful 
representation of the state of redshift information. 


C. LUMINOSITY FUNCTION AS TABLE0 

Available online at http://www.mpe.mpg.de/ 
jbuchner/agn_torus/paper/fieldsummary.cds 


9 SUPPLEMENTAL ONLINE MATERIAL 








